Trabajo Práctico Final: Análisis de Producción Lechera

Cátedra: Ciencia de Datos para la Toma de Decisiones

Profesor: Carlos Eduardo Marcos

Alumna: Luciana Inés Frías

Modelo aplicado: Random Forest

Universidad Tecnológica Nacional – Facultad Regional Trenque Lauquen

Objetivo del trabajo

El objetivo de este trabajo es construir un modelo predictivo que, a partir de la información recolectada durante los primeros dos meses de lactancia de cada vaca, permita estimar si esa vaca tendrá una producción lechera alta o baja en los 305 días posteriores al parto.

Para lograrlo, el modelo se alimenta de un conjunto de variables disponibles de forma temprana en el ciclo productivo, como:

  • La cantidad promedio diaria de leche en los primeros dos meses,

  • El porcentaje de grasa y proteínas en la leche,

  • El recuento de células somáticas (RCS) como indicador de salud mamaria,

  • La variación entre el primer y segundo mes en esos parámetros (litros, grasa, proteína, RCS),

  • La cantidad de partos (número de lactancias previas),

  • La producción total de leche acumulada durante la vida productiva del animal,

  • Y la cantidad de días transcurridos desde el último parto hasta el momento del control.

El propósito del modelo es poder tomar decisiones anticipadas sobre el manejo individual de cada vaca, aún sin esperar a completar toda la lactancia. Esto permitirá:

  • Detectar a tiempo vacas con alto potencial productivo, para enfocar recursos en ellas,

  • Optimizar la eficiencia del rodeo, priorizando alimentación, sanidad y decisiones reproductivas,

  • Y también identificar vacas con bajo rendimiento esperado, para ajustar su manejo o planificar su reemplazo.

En definitiva, este modelo busca ser una herramienta predictiva útil para mejorar la gestión del rodeo lechero desde una etapa muy temprana, basándose en datos reales, fácilmente disponibles, y orientados a potenciar la producción y rentabilidad del sistema.

1. Librerias

library(readr)
library(dplyr)
library(stringr)
library(knitr)
library(kableExtra)
library(tidyr)
library(ggplot2)
library(lubridate) 
library(readxl)
library(caret)
library(randomForest)

2. Carga de datos

En este caso utilizamos tres archivos

  1. datosproductivosbovinos: Contiene datos correspondiente a la fecha feb-mar 2025

  2. abril-mayo: Contiene datos correspondiente a la fecha abril-mayo 2024

  3. junio-julio: Contiene datos correspondiente a la fecha jun-jul 2024

Del primer archivo obtenemos la variable objetivo, correspondiente a “le_305_acu” la cual me define si la vaca fue de alta o de baja producción.

Del segundo y tercer archivo obtenemos las variables predictoras, que corresponden a datos obtenidos en los primeros dos meses posteriores al inicio de la lactancia de la vaca.

df <- read_delim("datosproductivosbovinos.csv", 
                                      delim = ";", escape_double = FALSE, locale = locale(decimal_mark = ",", 
                                                                                          grouping_mark = "."), trim_ws = TRUE)

abril_mayo <- read_excel("abril-mayo.xlsx")

junio_julio <- read_excel("junio-julio.xlsx")

3. Análisis del dataset primitivo

3.1. Estandarizar encabezados:

Colocar todo en minúscula y reemplazar todos los separadores (…/./-/espacio) a guión bajo”_”

colnames(df) <- tolower(colnames(df))
names(df) <- names(df) %>%
  str_replace_all("\\.\\.\\.", "_") %>%       # “...” → “_”
  str_replace_all("[\\.\\- ]", "_")           # “.”, “-” o espacio → “_”

3.2. Renombrar columnas para mejor lectura

El objetivo es colocar nombres más intuitivos para entender mejor de qué se tratan las variables “rcs_45…rcs_67”

df_modificado <- df %>%
  rename(
    rcs_abr_24  = `rcs_45`,
    rcs_may_24  = `rcs_47`,
    rcs_jun_24  = `rcs_49`,
    rcs_jul_24  = `rcs_51`,
    rcs_ago_24  = `rcs_53`,
    rcs_sept_24 = `rcs_55`,
    rcs_oct_24  = `rcs_57`,
    rcs_nov_24  = `rcs_59`,
    rcs_dic_24  = `rcs_61`,
    rcs_ene_25  = `rcs_63`,
    rcs_feb_25  = `rcs_65`,
    rcs_mar_25  = `rcs_67`,
    `sept_24`          = `set/24`
  )

3.3. Eliminar variables que no son de utilidad en el análisis

Variables que no aportan valor al análisis.

df_modificado <- df_modificado %>%
  select(
    -rp_madre,
    -nro_reg_madre,
    -f_servicio,
    -nro_servicio,
    -le_uc,
    -le_ca,
    -lote_ant,
    -so_uc,
    -so_ca,
    -re_ca,
    -dias_desde_ult_serv,
    -re_uc
  )

3.4. Diccionario de variables

Se construye el diccionario del dataset primitivo “datosbovinos” para entender mejor el significado de sus variables.

diccionario <- data.frame(
  variable_original = c(
    "estado", "f_nac", "padre", "est_actual", "nro_lact", "ini_lact",
    "secado", "apodo_toro", "prox_parto", "prox_secado",
    "gr_uc", "pr_uc", "gr_ca", "pr_ca",
    "dias_acu", "le_acu", "gr_acu", "pr_acu",
    "le_305_acu", "gr_305_acu", "pr_305_acu",
    "dias_vital", "le_vital", "gr_vital", "pr_vital",
    "lote_act", "tacto", "%_dif_uc/ca", "dias_al_parto",
    "f_fin_lact",
    "abr_24, may_24, jun_24,...,mar_25",
    "rcs_abr_24, rcs_may_24, rcs_jun_24,...,rcs_mar_25"
  ),
  nombre_variable = c(
    "estado fisiológico", "fecha de nacimiento", "padre de la vaca",
    "estado actual", "número de lactancias", "inicio de lactancia",
    "secado", "nombre del toro", "fecha de próximo parto",
    "fecha de secado", "% grasa marzo", "% proteína marzo",
    "% grasa febrero", "% proteína febrero", "días acumulados",
    "litros de leche acumulados", "kg de grasa acumulada",
    "kg de proteínas acumulada", "litros acumulados en 305 días",
    "kg grasa", "kg proteínas", "cantidad de días",
    "litros de leche", "kg de grasa", "kg de proteínas",
    "lote actual", "estado de preñez", "diferencia porcentual",
    "días acumulados desde el parto", "fin de lactancia",
    "litros de leche", "recuento de célula somática"
  ),
  tipo = c(
    "categórico (chr)", "fecha (Date)", "categórico (chr)",
    "categórico (chr)", "entero (int)", "fecha (Date)",
    "fecha (Date)", "categórico (chr)", "fecha (Date)",
    "fecha (Date)", "numérico (dbl)", "numérico (dbl)",
    "numérico (dbl)", "numérico (dbl)", "numérico (dbl)",
    "numérico (dbl)", "numérico (dbl)", "numérico (dbl)",
    "numérico (dbl)", "numérico (dbl)", "numérico (dbl)",
    "numérico (dbl)", "numérico (dbl)", "numérico (dbl)",
    "numérico (dbl)", "numérico (dbl)", "categórico (chr)",
    "numérico (dbl)", "numérico (dbl)", "fecha (Date)",
    "numérico (dbl)", "numérico (dbl)"
  ),
  descripción = c(
    "estado fisiológico general de la vaca (A = activa / R = rechazada).",
    "fecha de nacimiento del animal.",
    "identificador o nombre del toro padre.",
    "condición productiva al momento del registro (O = ordeñe / S = secado).",
    "cantidad de lactancias en curso (paridad).",
    "fecha del inicio de la lactancia del último parto.",
    "fecha en que se retira la vaca del ordeñe (≈60 días antes del próximo parto).",
    "identificador o nombre del toro utilizado para preñar a la vaca.",
    "fecha estimada del próximo parto.",
    "fecha programada de secado (próximo parto − 60 días).",
    "% de grasa butírica medido en marzo.",
    "% de proteínas medido en marzo.",
    "% de grasa butírica medido en febrero.",
    "% de proteínas medido en febrero.",
    "días acumulados de la lactancia actual hasta el control (07/03/2025).",
    "litros de leche acumulados en la lactancia actual.",
    "kilogramos de grasa acumulada en la lactancia actual.",
    "kilogramos de proteínas acumulada en la lactancia actual.",
    "litros acumulados en los 305 días posteriores al inicio de lactancia.",
    "kilogramos de grasa acumulados en los 305 días posteriores al inicio.",
    "kilogramos de proteínas acumulados en los 305 días posteriores al inicio.",
    "días totales de lactancia desde el primer parto.",
    "litros acumulados en su vida productiva.",
    "kilogramos de grasa producidos durante toda su vida productiva.",
    "kilogramos de proteínas producidos durante toda su vida productiva.",
    "lote o grupo productivo donde se encuentra la vaca.",
    "resultado del tacto rectal (P = preñada / V = vacía).",
    "diferencia porcentual de litros entre marzo y febrero.",
    "días transcurridos desde el último parto al registro (26/03/2025).",
    "fecha real de fin de lactancia (retiro del ordeñe).",
    "producción mensual de leche (L) por animal para cada mes indicado.",
    "valor de recuento de células somáticas mensual (células/ml)."
  ),
  stringsAsFactors = FALSE
)

A continuación se mostrará el diccionario realizado en formato tabla.

kable(diccionario, "html", escape = FALSE, align = "l") %>%
  kable_styling(
    full_width = FALSE,
    bootstrap_options = c("striped", "hover", "condensed")
  )
variable_original nombre_variable tipo descripción
estado estado fisiológico categórico (chr) estado fisiológico general de la vaca (A = activa / R = rechazada).
f_nac fecha de nacimiento fecha (Date) fecha de nacimiento del animal.
padre padre de la vaca categórico (chr) identificador o nombre del toro padre.
est_actual estado actual categórico (chr) condición productiva al momento del registro (O = ordeñe / S = secado).
nro_lact número de lactancias entero (int) cantidad de lactancias en curso (paridad).
ini_lact inicio de lactancia fecha (Date) fecha del inicio de la lactancia del último parto.
secado secado fecha (Date) fecha en que se retira la vaca del ordeñe (≈60 días antes del próximo parto).
apodo_toro nombre del toro categórico (chr) identificador o nombre del toro utilizado para preñar a la vaca.
prox_parto fecha de próximo parto fecha (Date) fecha estimada del próximo parto.
prox_secado fecha de secado fecha (Date) fecha programada de secado (próximo parto − 60 días).
gr_uc % grasa marzo numérico (dbl) % de grasa butírica medido en marzo.
pr_uc % proteína marzo numérico (dbl) % de proteínas medido en marzo.
gr_ca % grasa febrero numérico (dbl) % de grasa butírica medido en febrero.
pr_ca % proteína febrero numérico (dbl) % de proteínas medido en febrero.
dias_acu días acumulados numérico (dbl) días acumulados de la lactancia actual hasta el control (07/03/2025).
le_acu litros de leche acumulados numérico (dbl) litros de leche acumulados en la lactancia actual.
gr_acu kg de grasa acumulada numérico (dbl) kilogramos de grasa acumulada en la lactancia actual.
pr_acu kg de proteínas acumulada numérico (dbl) kilogramos de proteínas acumulada en la lactancia actual.
le_305_acu litros acumulados en 305 días numérico (dbl) litros acumulados en los 305 días posteriores al inicio de lactancia.
gr_305_acu kg grasa numérico (dbl) kilogramos de grasa acumulados en los 305 días posteriores al inicio.
pr_305_acu kg proteínas numérico (dbl) kilogramos de proteínas acumulados en los 305 días posteriores al inicio.
dias_vital cantidad de días numérico (dbl) días totales de lactancia desde el primer parto.
le_vital litros de leche numérico (dbl) litros acumulados en su vida productiva.
gr_vital kg de grasa numérico (dbl) kilogramos de grasa producidos durante toda su vida productiva.
pr_vital kg de proteínas numérico (dbl) kilogramos de proteínas producidos durante toda su vida productiva.
lote_act lote actual numérico (dbl) lote o grupo productivo donde se encuentra la vaca.
tacto estado de preñez categórico (chr) resultado del tacto rectal (P = preñada / V = vacía).
%_dif_uc/ca diferencia porcentual numérico (dbl) diferencia porcentual de litros entre marzo y febrero.
dias_al_parto días acumulados desde el parto numérico (dbl) días transcurridos desde el último parto al registro (26/03/2025).
f_fin_lact fin de lactancia fecha (Date) fecha real de fin de lactancia (retiro del ordeñe).
abr_24, may_24, jun_24,…,mar_25 litros de leche numérico (dbl) producción mensual de leche (L) por animal para cada mes indicado.
rcs_abr_24, rcs_may_24, rcs_jun_24,…,rcs_mar_25 recuento de célula somática numérico (dbl) valor de recuento de células somáticas mensual (células/ml).

4. Exporación inicial de datos (EDA)

4.1. Analizar valores faltantes:

colSums(is.na(df_modificado))
##            rp        estado         f_nac         padre    est_actual 
##             1             0             0             8             0 
##      nro_lact      ini_lact        secado    apodo_toro    prox_parto 
##             0             0             0           262             0 
##   prox_secado         gr_uc         pr_uc         gr_ca         pr_ca 
##             0           298           298           369           369 
##      dias_acu        le_acu        gr_acu        pr_acu    le_305_acu 
##             0            47            47            47           877 
##    gr_305_acu    pr_305_acu    dias_vital      le_vital      gr_vital 
##           877           877             0             0             0 
##      pr_vital      lote_act         tacto   %_dif_uc/ca dias_al_parto 
##             0             0           475             0             0 
##    f_fin_lact        abr_24    rcs_abr_24        may_24    rcs_may_24 
##             0           496           496           462           462 
##        jun_24    rcs_jun_24        jul_24    rcs_jul_24        ago_24 
##           354           354           381           381           245 
##    rcs_ago_24       sept_24   rcs_sept_24        oct_24    rcs_oct_24 
##           245           195           195           107           107 
##        nov_24    rcs_nov_24        dic_24    rcs_dic_24        ene_25 
##            74            74            55            55           183 
##    rcs_ene_25        feb_25    rcs_feb_25        mar_25    rcs_mar_25 
##           183           248           248           214           214

La inspección de los registros revela que la mayoría de los valores ausentes se explican por cuestiones temporales y no por errores de captura:

  1. grasa, proteinas marzo=298 (esto corresponde a que la vaca no estuvo en ordeñe en ese mes, es por esto que aparece faltante el dato)

  2. grasa, proteinas feb=369 (esto corresponde a que la vaca no estuvo en ordeñe en ese mes, es por esto que aparece faltante el dato)

  3. Le.ACU,Gr.ACU, Pr.ACU=47 (la vaca inicia la lactancia posterior a la fecha decontrol del tambero [07/03/2025])

  4. Le.305 ACU,Gr.305 ACU, Pr.305 ACU=877 (esto indica que la vaca aún no superó los 305 días posteriores a su lactancia, por eso no aparece el dato)

  5. mes-año y rcs-mes-año (el valor faltante indica que la lactancia comenzó posterior al mes-año señalado, por eso no está el dato)

  6. padre=8 (se tratará como desconocido)

  7. apodo toro=262 (son vacas que no han sido inseminadas para preñez). serán tradas como NC (no corresp.)

  8. tacto=475 (vacas inseminadas recientemente, en teoria debe superar entre los 28 y 45 días post-inseminación) tratadas como NC

4.2. Reemplazo de los valores faltantes en variables categóricas

En el caso de tacto y apodo toro, los valores faltantes pertenecen a animales aún no han sido inseminados o están en ventana post-inseminación; se etiquetan como “no aplica”.

vars_reemplazo <- c("padre", "apodo_toro", "tacto")  

df_modificado <- df_modificado %>%
  mutate(across(
    all_of(vars_reemplazo),
    ~ replace_na(.x, "no aplica")
  ))

Estructura y resumen estadístico

glimpse(df_modificado)   
## Rows: 1,318
## Columns: 55
## $ rp            <dbl> 1, 4, 5, 18, 19, 23, 29, 30, 32, 36, 40, 45, 55, 56, 64,…
## $ estado        <chr> "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "…
## $ f_nac         <chr> "5/11/2022", "12/11/2022", "8/11/2022", "18/11/2022", "1…
## $ padre         <chr> "DELTA-F", "DELTA-F", "DAVOS", "FACTUAL", "Revelat", "Re…
## $ est_actual    <chr> "O", "O", "O", "O", "O", "O", "O", "O", "O", "O", "O", "…
## $ nro_lact      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ ini_lact      <chr> "4/9/2024", "19/8/2024", "31/8/2024", "12/9/2024", "6/9/…
## $ secado        <chr> "/  /", "/  /", "/  /", "/  /", "/  /", "/  /", "/  /", …
## $ apodo_toro    <chr> "FOREVER", "Kareem", "CATFISH", "no aplica", "FOREVER", …
## $ prox_parto    <chr> "/  /", "3/9/2025", "/  /", "/  /", "/  /", "/  /", "/  …
## $ prox_secado   <chr> "/  /", "5/7/2025", "/  /", "/  /", "/  /", "/  /", "/  …
## $ gr_uc         <dbl> 2.60, 3.38, 4.27, 1.50, 2.77, 2.56, 4.01, 2.87, 3.28, 3.…
## $ pr_uc         <dbl> 3.21, 3.04, 3.34, 3.27, 4.00, 3.24, 3.59, 3.47, 3.32, 3.…
## $ gr_ca         <dbl> 3.38, 2.98, 3.55, 3.10, 3.45, 3.47, 3.78, 3.47, 2.50, 4.…
## $ pr_ca         <dbl> 3.10, 3.02, 3.52, 3.05, 3.21, 3.56, 3.56, 3.35, 3.40, 3.…
## $ dias_acu      <dbl> 185, 201, 189, 177, 183, 47, 186, 185, 178, 182, 174, 18…
## $ le_acu        <dbl> 6713.0, 7611.3, 6600.3, 7874.5, 8011.1, 1543.2, 7033.2, …
## $ gr_acu        <dbl> 236.34, 208.96, 209.83, 199.53, 206.44, 48.35, 228.74, 2…
## $ pr_acu        <dbl> 215.46, 262.44, 247.49, 247.25, 274.95, 53.17, 251.12, 2…
## $ le_305_acu    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ gr_305_acu    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ pr_305_acu    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ dias_vital    <dbl> 185, 201, 189, 177, 183, 47, 186, 185, 178, 182, 174, 18…
## $ le_vital      <dbl> 6713.0, 7611.3, 6600.3, 7874.5, 8011.1, 1543.2, 7033.2, …
## $ gr_vital      <dbl> 236.34, 208.96, 209.83, 199.53, 206.44, 48.35, 228.74, 2…
## $ pr_vital      <dbl> 215.46, 262.44, 247.49, 247.25, 274.95, 53.17, 251.12, 2…
## $ lote_act      <dbl> 7, 5, 6, 5, 3, 7, 7, 6, 3, 3, 3, 7, 2, 6, 6, 7, 6, 3, 6,…
## $ tacto         <chr> "no aplica", "P", "V", "V", "no aplica", "no aplica", "n…
## $ `%_dif_uc/ca` <chr> "2.73", "-1.49", "-20.21", "0", "2.44", "22.88", "16.07"…
## $ dias_al_parto <dbl> 203, 219, 207, 195, 201, 65, 204, 203, 196, 200, 192, 19…
## $ f_fin_lact    <chr> "/  /", "/  /", "/  /", "/  /", "/  /", "/  /", "/  /", …
## $ abr_24        <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ rcs_abr_24    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ may_24        <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ rcs_may_24    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ jun_24        <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ rcs_jun_24    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ jul_24        <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ rcs_jul_24    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ ago_24        <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ rcs_ago_24    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ sept_24       <dbl> NA, 26.2, 20.0, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ rcs_sept_24   <dbl> NA, 403, 403, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ oct_24        <dbl> 31.6, 36.0, 31.2, 37.2, 40.2, NA, 34.8, 38.4, 40.0, 35.4…
## $ rcs_oct_24    <dbl> 2162, 34, 59, 94, 10, NA, 14, 16, 30, 19, 21, 16, 41, 9,…
## $ nov_24        <dbl> 35.0, 40.0, 34.6, 50.6, 43.6, NA, 40.0, 41.0, 51.8, 40.4…
## $ rcs_nov_24    <dbl> 603, 22, 76, 0, 30, NA, 9, 4, 323, 8, 4, 27, 376, 20, 24…
## $ dic_24        <dbl> 41.0, 43.6, 43.0, 48.0, 53.0, NA, 41.6, 47.0, 61.6, 52.0…
## $ rcs_dic_24    <dbl> 355, 6, 105, 11, 7, NA, 7, 14, 20, 15, 6, 23, 64, 12, 0,…
## $ ene_25        <dbl> 38.6, 40.4, 39.6, 48.0, 43.0, NA, 39.6, 44.6, 49.2, 45.4…
## $ rcs_ene_25    <dbl> 172, 12, 184, 88, 8, NA, 5, 5, 7, 6, 15, 25, 41, 5, 7, 6…
## $ feb_25        <dbl> 36.6, 40.2, 37.6, 41.4, 41.0, 30.6, 33.6, 41.0, 46.6, 43…
## $ rcs_feb_25    <dbl> 560, 174, 85, 9, 2, 236, 12, 19, 19, 7, 7, 19, 51, 3, 4,…
## $ mar_25        <dbl> 37.6, 39.6, 30.0, 41.4, 42.0, 37.6, 39.0, 42.4, 49.6, 45…
## $ rcs_mar_25    <dbl> 74, 9, 63, 5, 5, 11, 6, 12, 240, 12, 16, 165, 5, 9, 16, …
summary(df_modificado)
##        rp           estado             f_nac              padre          
##  Min.   :    1   Length:1318        Length:1318        Length:1318       
##  1st Qu.:16570   Class :character   Class :character   Class :character  
##  Median :17805   Mode  :character   Mode  :character   Mode  :character  
##  Mean   :17349                                                           
##  3rd Qu.:18821                                                           
##  Max.   :19995                                                           
##  NA's   :1                                                               
##   est_actual           nro_lact       ini_lact            secado         
##  Length:1318        Min.   :1.000   Length:1318        Length:1318       
##  Class :character   1st Qu.:1.000   Class :character   Class :character  
##  Mode  :character   Median :2.000   Mode  :character   Mode  :character  
##                     Mean   :1.869                                        
##                     3rd Qu.:3.000                                        
##                     Max.   :5.000                                        
##                                                                          
##   apodo_toro         prox_parto        prox_secado            gr_uc      
##  Length:1318        Length:1318        Length:1318        Min.   :1.130  
##  Class :character   Class :character   Class :character   1st Qu.:2.558  
##  Mode  :character   Mode  :character   Mode  :character   Median :3.290  
##                                                           Mean   :3.328  
##                                                           3rd Qu.:3.960  
##                                                           Max.   :9.520  
##                                                           NA's   :298    
##      pr_uc           gr_ca           pr_ca          dias_acu    
##  Min.   :2.670   Min.   :1.040   Min.   :2.580   Min.   :  0.0  
##  1st Qu.:3.270   1st Qu.:2.800   1st Qu.:3.220   1st Qu.:181.0  
##  Median :3.450   Median :3.290   Median :3.400   Median :260.0  
##  Mean   :3.537   Mean   :3.252   Mean   :3.435   Mean   :251.1  
##  3rd Qu.:3.670   3rd Qu.:3.740   3rd Qu.:3.600   3rd Qu.:334.0  
##  Max.   :5.740   Max.   :7.630   Max.   :5.470   Max.   :696.0  
##  NA's   :298     NA's   :369     NA's   :369                    
##      le_acu          gr_acu            pr_acu          le_305_acu   
##  Min.   :   48   Min.   :   1.66   Min.   :   2.28   Min.   : 6720  
##  1st Qu.: 8108   1st Qu.: 225.41   1st Qu.: 286.70   1st Qu.:12899  
##  Median :11750   Median : 322.59   Median : 410.58   Median :14113  
##  Mean   :11564   Mean   : 337.35   Mean   : 405.23   Mean   :14060  
##  3rd Qu.:15336   3rd Qu.: 438.37   3rd Qu.: 541.83   3rd Qu.:15710  
##  Max.   :32747   Max.   :1037.32   Max.   :1108.77   Max.   :18698  
##  NA's   :47      NA's   :47        NA's   :47        NA's   :877    
##    gr_305_acu      pr_305_acu      dias_vital        le_vital    
##  Min.   :239.2   Min.   :248.4   Min.   :   5.0   Min.   :   70  
##  1st Qu.:362.9   1st Qu.:448.4   1st Qu.: 296.0   1st Qu.:11430  
##  Median :421.4   Median :497.8   Median : 486.0   Median :19998  
##  Mean   :424.5   Mean   :491.6   Mean   : 548.2   Mean   :23536  
##  3rd Qu.:481.6   3rd Qu.:543.9   3rd Qu.: 796.0   3rd Qu.:34041  
##  Max.   :687.0   Max.   :642.0   Max.   :1547.0   Max.   :80670  
##  NA's   :877     NA's   :877                                     
##     gr_vital          pr_vital          lote_act        tacto          
##  Min.   :   3.21   Min.   :   2.89   Min.   :1.000   Length:1318       
##  1st Qu.: 340.46   1st Qu.: 409.44   1st Qu.:3.000   Class :character  
##  Median : 667.29   Median : 709.14   Median :4.000   Mode  :character  
##  Mean   : 778.79   Mean   : 820.39   Mean   :4.268                     
##  3rd Qu.:1158.16   3rd Qu.:1192.50   3rd Qu.:6.000                     
##  Max.   :2780.31   Max.   :2718.52   Max.   :7.000                     
##                                                                        
##  %_dif_uc/ca        dias_al_parto    f_fin_lact            abr_24     
##  Length:1318        Min.   :  8.0   Length:1318        Min.   : 4.00  
##  Class :character   1st Qu.:199.0   Class :character   1st Qu.:37.00  
##  Mode  :character   Median :281.0   Mode  :character   Median :41.30  
##                     Mean   :275.1                      Mean   :40.87  
##                     3rd Qu.:367.0                      3rd Qu.:46.00  
##                     Max.   :754.0                      Max.   :67.00  
##                                                        NA's   :496    
##    rcs_abr_24         may_24        rcs_may_24         jun_24     
##  Min.   :   0.0   Min.   : 2.00   Min.   :   0.0   Min.   : 2.00  
##  1st Qu.:  15.0   1st Qu.:39.00   1st Qu.:   8.0   1st Qu.:37.00  
##  Median :  40.0   Median :45.00   Median :  15.0   Median :43.00  
##  Mean   : 178.9   Mean   :43.32   Mean   : 165.2   Mean   :42.28  
##  3rd Qu.: 145.8   3rd Qu.:51.00   3rd Qu.: 129.0   3rd Qu.:49.65  
##  Max.   :5188.0   Max.   :76.00   Max.   :5393.0   Max.   :72.80  
##  NA's   :496      NA's   :462     NA's   :462      NA's   :354    
##    rcs_jun_24          jul_24        rcs_jul_24         ago_24     
##  Min.   :   0.00   Min.   :15.00   Min.   :   0.0   Min.   : 3.00  
##  1st Qu.:   7.75   1st Qu.:42.20   1st Qu.:   8.0   1st Qu.:41.00  
##  Median :  14.00   Median :48.00   Median :  17.0   Median :47.00  
##  Mean   : 173.53   Mean   :48.03   Mean   : 132.3   Mean   :46.68  
##  3rd Qu.: 123.25   3rd Qu.:56.60   3rd Qu.: 102.0   3rd Qu.:55.00  
##  Max.   :5159.00   Max.   :87.20   Max.   :4141.0   Max.   :82.00  
##  NA's   :354       NA's   :381     NA's   :381      NA's   :245    
##    rcs_ago_24        sept_24       rcs_sept_24         oct_24     
##  Min.   :   0.0   Min.   :11.00   Min.   :   0.0   Min.   : 8.00  
##  1st Qu.:   9.0   1st Qu.:40.60   1st Qu.:   6.0   1st Qu.:38.00  
##  Median :  19.0   Median :46.00   Median :  16.0   Median :43.00  
##  Mean   : 107.2   Mean   :45.37   Mean   : 153.5   Mean   :43.02  
##  3rd Qu.:  98.0   3rd Qu.:52.80   3rd Qu.: 118.5   3rd Qu.:49.40  
##  Max.   :4599.0   Max.   :82.00   Max.   :3877.0   Max.   :75.60  
##  NA's   :245      NA's   :195     NA's   :195      NA's   :107    
##    rcs_oct_24         nov_24        rcs_nov_24         dic_24     
##  Min.   :   0.0   Min.   : 5.00   Min.   :   0.0   Min.   : 3.00  
##  1st Qu.:   9.0   1st Qu.:37.60   1st Qu.:  10.0   1st Qu.:41.60  
##  Median :  26.0   Median :43.50   Median :  23.0   Median :47.40  
##  Mean   : 110.9   Mean   :41.07   Mean   : 111.3   Mean   :45.25  
##  3rd Qu.:  75.5   3rd Qu.:49.40   3rd Qu.:  95.0   3rd Qu.:53.60  
##  Max.   :3713.0   Max.   :70.60   Max.   :3369.0   Max.   :80.60  
##  NA's   :107      NA's   :74      NA's   :74       NA's   :55     
##    rcs_dic_24         ene_25        rcs_ene_25         feb_25     
##  Min.   :   0.0   Min.   : 4.00   Min.   :   0.0   Min.   : 4.00  
##  1st Qu.:  10.0   1st Qu.:37.00   1st Qu.:   6.0   1st Qu.:34.00  
##  Median :  37.0   Median :42.60   Median :  15.0   Median :39.60  
##  Mean   : 141.7   Mean   :42.13   Mean   : 100.2   Mean   :38.63  
##  3rd Qu.: 106.5   3rd Qu.:49.20   3rd Qu.:  67.0   3rd Qu.:44.40  
##  Max.   :4740.0   Max.   :69.60   Max.   :4315.0   Max.   :62.00  
##  NA's   :55       NA's   :183     NA's   :183      NA's   :248    
##    rcs_feb_25         mar_25        rcs_mar_25    
##  Min.   :   0.0   Min.   : 2.00   Min.   :   0.0  
##  1st Qu.:   8.0   1st Qu.:31.40   1st Qu.:   9.0  
##  Median :  29.0   Median :39.40   Median :  41.5  
##  Mean   : 158.4   Mean   :36.73   Mean   : 181.6  
##  3rd Qu.: 102.8   3rd Qu.:44.00   3rd Qu.: 136.0  
##  Max.   :5551.0   Max.   :58.20   Max.   :4126.0  
##  NA's   :248      NA's   :214     NA's   :214

4.3. Análisis univariado

1. Leche acumulada desde el inicio de lactancia

#Histograma
ggplot(df_modificado, aes(x = le_acu)) +
  geom_histogram(bins = 30, fill = "#4DA8DA", colour = "white") +
  labs(title = "Cantidad de litros acum en lact. actual",
       x = "Litros acumulados", y = "Frecuencia") +
  theme_minimal()

# Boxplot
ggplot(df_modificado, aes(y = le_acu)) +
  geom_boxplot(fill = "#FFBC42", colour = "black") +
  labs(title = "Boxplot de litros acum en lact. actual",
       y = "Litros acumulados") +
  theme_minimal()

La producción acumulada (Le.ACU) muestra una distribución asimétrica a la derecha, con la mayoría de las vacas entre 7 000 y 18 000 L y algunos casos extremos que llegan a 32 000 L, asociados a lactancias prolongadas. La mediana es de ~11 500 L, con outliers por encima de 25 000 L que se deben a vacas que tuvieron abortos y por eso siguen lactando, pero no se registra la fecha de tal aborto.

2. Kg de grasa acumulada desde el inicio de la lactancia

# Histograma
ggplot(df_modificado, aes(x = gr_acu)) +
  geom_histogram(bins = 30, fill = "#4DA8DA", colour = "white") +
  labs(title = "Distribución de Gr.ACU",
       x = "Kg de grasa acumulada", y = "Frecuencia") +
  theme_minimal()

# Boxplot
ggplot(df_modificado, aes(y = gr_acu)) +
  geom_boxplot(fill = "#FFBC42", colour = "black") +
  labs(title = "Boxplot de Gr.ACU",
       y = "Kg de grasa acumulada") +
  theme_minimal()

La producción de grasa acumulada muestra una distribución unimodal con pico en 300 kg y mayoría entre 200 y 450 kg. La mediana es de ~330 kg, con outliers que superan los 750 kg, asociados a lactancias largas o persistentes.

3. Kg de proteína acumulada

# Histograma
ggplot(df_modificado, aes(x = pr_acu)) +
  geom_histogram(bins = 30, fill = "#4DA8DA", colour = "white") +
  labs(title = "Distribución de Pr.ACU",
       x = "Kg de proteína acumulada", y = "Frecuencia") +
  theme_minimal()

# Boxplot
ggplot(df_modificado, aes(y = pr_acu)) +
  geom_boxplot(fill = "#FFBC42", colour = "black") +
  labs(title = "Boxplot de Pr.ACU",
       y = "Kg de proteína acumulada") +
  theme_minimal()

La producción de proteína acumulada es asimétrica a la derecha, con mayoría entre 250 y 550 kg y un pico en 350 kg. La mediana es ~370 kg, y los outliers superan los 750 kg, vinculados a vacas con lactancias largas o muy productivas. El rodeo es homogéneo, salvo por un grupo reducido de alto rendimiento.

4. Cantidad de litros acumulados en los primeros 305 días posterior a la lactancia

# Histograma
ggplot(df_modificado, aes(x = `le_305_acu`)) +
  geom_histogram(bins = 30, fill = "#4DA8DA", colour = "white") +
  labs(title = "Distribución de Le.305 ACU",
       x = "Litros proyectados a 305 días", y = "Frecuencia") +
  theme_minimal()

# Boxplot
ggplot(df_modificado, aes(y = `le_305_acu`)) +
  geom_boxplot(fill = "#FFBC42", colour = "black") +
  labs(title = "Boxplot de Le.305 ACU",
       y = "Litros proyectados a 305 días") +
  theme_minimal()

La producción proyectada a 305 días se concentra entre 13 000 y 16 000 L, con una mediana de ~14 500 L y variabilidad moderada. Hay pocos outliers bajos (menos de 9 000 L), asociados a vacas con lactancias recientes o datos incompletos, y casi no se registran valores extremos altos, lo que indica un potencial productivo bien delimitado.

5. Cantidad de litros producidos desde su primer lactancia

# Histograma
ggplot(df_modificado, aes(x = le_vital)) +
  geom_histogram(bins = 30, fill = "#4DA8DA", colour = "white") +
  labs(title = "Distribución de Le.Vital",
       x = "Litros vitalicios", y = "Frecuencia") +
  theme_minimal()

# Boxplot
ggplot(df_modificado, aes(y = le_vital)) +
  geom_boxplot(fill = "#FFBC42", colour = "black") +
  labs(title = "Boxplot de Le.Vital",
       y = "Litros vitalicios") +
  theme_minimal()

La producción vitalicia muestra una distribución sesgada a la derecha, con mayoría entre 10 000 y 35 000 L y un pico entre 15 000 y 20 000 L. La mediana es de ~20 000 L y el rango es amplio, reflejando alta variabilidad. Existen outliers por encima de 75 000 L, correspondientes a vacas con vidas productivas muy largas o datos extremos, y algunos valores bajos (casi 0) por lactancias recientes o registros incompletos.

6. % de grasa en el mes de marzo

# Histograma
ggplot(df_modificado, aes(x = `gr_uc`)) +
  geom_histogram(bins = 30, fill = "#4DA8DA", colour = "white") +
  labs(title = "Distribución de grasa-marzo",
       x = "% Grasa en marzo", y = "Frecuencia") +
  theme_minimal()

# Boxplot
ggplot(df_modificado, aes(y = `gr_uc`)) +
  geom_boxplot(fill = "#FFBC42", colour = "black") +
  labs(title = "Boxplot de grasa-marzo",
       y = "% Grasa en marzo") +
  theme_minimal()

La grasa en leche de marzo presenta una distribución casi simétrica, con la mayoría de las vacas entre 2.5 % y 4.5 % y un pico en 3 %–3.5 %. La mediana es de ~3.4 %, con variabilidad moderada. Los outliers por encima de 5.5 % podrían corresponder a vacas con alta composición grasa o posibles errores. Valores bajos (<2 %) son poco frecuentes.

7. % proteínas en el mes de febrero

# Histograma
ggplot(df_modificado, aes(x = `pr_uc`)) +
  geom_histogram(bins = 30, fill = "#4DA8DA", colour = "white") +
  labs(title = "Distribución de proteinas-marzo",
       x = "% Proteína en marzo", y = "Frecuencia") +
  theme_minimal()

# Boxplot
ggplot(df_modificado, aes(y = `pr_uc`)) +
  geom_boxplot(fill = "#FFBC42", colour = "black") +
  labs(title = "Boxplot de proteinas-marzo",
       y = "% Proteína en marzo") +
  theme_minimal()

La proteína en leche se concentra entre 3.2 % y 3.8 %, con un pico en 3.5 % y mediana de ~3.4 %, lo que refleja consistencia proteica en el rodeo. Hay pocos valores bajos (<3.0 %) y una cola derecha moderada con outliers por encima de 4.2 %, que podrían indicar composición inusualmente alta o posibles errores de laboratorio/genética especial.

8. Recuento de célula somáticas en los primeros tres meses de 2025

# Pivot a formato largo, filtrando solo meses 25
rcs_2025 <- df_modificado %>%
  select(matches("^rcs_.*_25$")) %>%   # rcs_ene_25, rcs_feb_25, rcs_mar_25
  pivot_longer(
    cols = everything(),
    names_to  = "mes",
    values_to = "rcs"
  )

# Histograma facetado
ggplot(rcs_2025, aes(x = rcs)) +
  geom_histogram(bins = 30, fill = "#4DA8DA", colour = "white") +
  facet_wrap(~ mes, scales = "free_x", ncol = 1) +
  labs(x = "RCS (células/ml)", y = "Frecuencia") +
  theme_minimal()

# Boxplot facetado
ggplot(rcs_2025, aes(x = mes, y = rcs)) +
  geom_boxplot(fill = "#FFBC42", colour = "black") +
  coord_flip() +
  labs(x = "", y = "RCS (células/ml)") +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 10))

El RCS en enero, febrero y marzo muestra que más del 60 % del rodeo tiene niveles muy bajos (<200 cél/ml), indicando buena salud general. Las distribuciones son similares, con colas largas hacia la derecha y algunos outliers elevados, especialmente en febrero, lo que sugiere un ligero aumento de mastitis ese mes. Las medianas están cerca de 0 y el rango intercuartílico es bajo, reflejando estabilidad en los valores.

9. Cantidad de días lactando

df_modificado %>%
  select(`dias_acu`, `dias_vital`) %>%
  pivot_longer(
    cols = everything(),
    names_to  = "tipo_dias",
    values_to = "dias"
  ) %>%
  ggplot(aes(x = tipo_dias, y = dias, fill = tipo_dias)) +
  geom_boxplot(colour = "black") +
  coord_flip() +
  labs(x = "", y = "Días") +
  theme_minimal() +
  theme(legend.position = "none")

Días ACU (lactancia acumulada): La mediana es de ~300 días, lo que indica que la mayoría de las vacas están próximas a completar el periodo ideal de lactancia de 305 días. Hay outliers que superan los 550 días, lo cual podría deberse a lactancias muy prolongadas o errores de registro.

Días Vital (edad productiva desde primer parto): Mediana de ~500 días, con gran variabilidad entre vacas jóvenes y veteranas. Algunos casos extremos superan los 1.000 días, posiblemente por múltiples lactancias o posibles abortos no registrados.

10. Producción mensual 2025

prod_2025 <- df_modificado %>%
  select(`ene_25`, `feb_25`, `mar_25`) %>%
  pivot_longer(
    cols = everything(),
    names_to  = "mes",
    values_to = "litros"
  )

# 2. Histograma facetado
ggplot(prod_2025, aes(x = litros)) +
  geom_histogram(bins = 30, fill = "#4DA8DA", colour = "white") +
  facet_wrap(~ mes, scales = "free_x", ncol = 1) +
  labs(x = "Litros por mes", y = "Frecuencia") +
  theme_minimal()

# 3. Boxplot facetado
ggplot(prod_2025, aes(x = mes, y = litros, fill = mes)) +
  geom_boxplot(colour = "black") +
  coord_flip() +
  labs(x = "", y = "Litros por mes") +
  theme_minimal() +
  theme(legend.position = "none")

Enero tuvo la mayor mediana (42.6 L) y más vacas de alta producción. En febrero y marzo la mediana bajó (~39.5 L) y hubo menos casos extremos, indicando un leve descenso estacional.

11. Estado productivo

# 1. Preparar datos en formato largo
estados_long <- df_modificado %>%
  select(estado, `est_actual`) %>%
  pivot_longer(
    cols = everything(),
    names_to  = "tipo_estado",
    values_to = "valor"
  )

# 2. Gráfico de barras facetado
ggplot(estados_long, aes(x = valor, fill = valor)) +
  geom_bar() +
  facet_wrap(~ tipo_estado, scales = "free_x", ncol = 1) +
  labs(x = "Categoría", y = "Cantidad de vacas") +
  theme_minimal() +
  theme(legend.position = "none")

La mayoría de las vacas están activas y en ordeñe. Para modelar producción, se recomienda enfocarse solo en vacas activas y en ordeñe, excluyendo las retiradas y secas.

12. Tacto (reproducción)

ggplot(df_modificado, aes(x = tacto, fill = tacto)) +
  geom_bar() +
  scale_fill_manual(values = c("P" = "#4DA8DA", "V" = "#FFBC42")) +
  labs(title = "Frecuencia de tacto: preñada vs vacía",
       x = "Resultado de tacto", y = "Cantidad de vacas") +
  theme_minimal() +
  theme(legend.position = "none")

La mayoría de las vacas están preñadas (~700), unas 150 están vacías y ~500 no aplican por estar secas o sin servicio reciente. Para evaluar eficiencia reproductiva, conviene excluir estas últimas.

13. Categorías de manejo

manejo_long <- df_modificado %>%
  select(`lote_act`, `nro_lact`) %>%
  mutate(
    `lote_act`   = as.character(`lote_act`),
    `nro_lact` = as.character(`nro_lact`)
  ) %>%
  pivot_longer(
    cols      = everything(),
    names_to  = "categoria",
    values_to = "valor"
  )

ggplot(manejo_long, aes(x = valor, fill = valor)) +
  geom_bar() +
  facet_wrap(~ categoria, scales = "free_x", ncol = 1) +
  labs(x = "Nivel", y = "Cantidad de vacas") +
  theme_minimal() +
  theme(legend.position = "none",
        strip.text = element_text(face = "bold"))

El lote 4 es el más grande y podría influir en los análisis. La mayoría de las vacas está en su primera o segunda lactancia, con pocas vacas veteranas.

4.4. Análisis bivariado

1. Le.ACU vs. Dias ACU: curva lactacional

ggplot(df_modificado, aes(x = `dias_acu`, y = le_acu)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", se = TRUE, colour = "#E15759") +
  labs(
    title = "Le.ACU vs. Dias ACU",
    x = "Días en lactancia (Dias ACU)",
    y = "Litros acumulados (Le.ACU)"
  ) +
  theme_minimal()

La producción acumulada sube rápido hasta ~200 días, se estabiliza entre 100–300 días y luego crece más lento. Outliers indican posibles errores de registro o problemas de salud/manejo.

2. Le.ACU vs. nro lactancia: producción por paridad

ggplot(df_modificado, aes(x = factor(`nro_lact`), y = le_acu, fill = factor(`nro_lact`))) +
  geom_boxplot() +
  geom_jitter(width = 0.2, alpha = 0.3, size = 0.8) +
  labs(
    title = "Le.ACU vs. Número de lactancia",
    x = "Paridad (nro lactancia)",
    y = "Litros acumulados (Le.ACU)"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

La mayor producción media se da en la 2ª y 5ª lactancia, aunque esta última está influida por pocas vacas superproductoras. La 1ª y 3ª tienen menor rendimiento y la 4ª muestra alta variabilidad.

3. Le.ACU vs. lote actual: producción por lote

ggplot(df_modificado, aes(x = factor(`lote_act`), y = le_acu, fill = factor(`lote_act`))) +
  geom_boxplot() +
  geom_jitter(width = 0.2, alpha = 0.3, size = 0.8) +
  labs(
    title = "Le.ACU vs. Lote actual",
    x = "Lote actual",
    y = "Litros acumulados (Le.ACU)"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

El lote 4 tiene la mayor producción y variabilidad. Los lotes 5, 6 y 7 le siguen con medianas similares (~12 000–13 000 L). Lote 1 tiene producción muy baja; probablemente incluye vacas no productivas.

4. cant de lts vs mes

# Vector con el orden deseado
meses <- c("abr_24","may_24","jun_24","jul_24","ago_24",
           "sept_24","oct_24","nov_24","dic_24",
           "ene_25","feb_25","mar_25")

# 1. Seleccionar columnas y convertir a numérico
df_promedios <- df_modificado %>%
  select(all_of(meses)) %>%
  mutate(across(everything(), ~ as.numeric(gsub(",", ".", .))))

# 2. Calcular medias y pivotar a formato largo
promedios <- df_promedios %>%
  summarise(across(everything(), ~ mean(.x, na.rm = TRUE))) %>%
  pivot_longer(
    everything(),
    names_to  = "Mes",
    values_to = "Promedio"
  ) %>%
  mutate(Mes = factor(Mes, levels = meses))  # asegurar orden
# 3. Mostrar tabla de promedios
print(promedios)
## # A tibble: 12 × 2
##    Mes     Promedio
##    <fct>      <dbl>
##  1 abr_24      40.9
##  2 may_24      43.3
##  3 jun_24      42.3
##  4 jul_24      48.0
##  5 ago_24      46.7
##  6 sept_24     45.4
##  7 oct_24      43.0
##  8 nov_24      41.1
##  9 dic_24      45.2
## 10 ene_25      42.1
## 11 feb_25      38.6
## 12 mar_25      36.7
# 4. Gráfico Mes vs. Litros promedio en orden cronológico
ggplot(promedios, aes(x = Mes, y = Promedio, group = 1)) +
  geom_line(colour = "#4DA8DA", size = 1) +
  geom_point(colour = "#E15759", size = 3) +
  labs(
    title = "Producción promedio diario (abr-24 a mar-25)",
    x = "Mes",
    y = "Litros promedio por día"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

La producción sube hasta julio, baja en primavera, repunta en diciembre y cae fuerte en verano. Máximo en invierno, mínimo en marzo por efecto estacional.

5. Preparación del dataset

Filtramos y dejamos únicamente los registros que tengan datos en le_305_acu, ya que esta variable me define si la vaca fue de alta o baja producción.

df_nuevo <- df_modificado[!is.na(df_modificado$le_305_acu), ]

Convertimos a fecha la variable, ya que también filtramos por vacas que iniciaron su lactancia posterior a marzo porque tenemos registros desde abril 2024 en adelante.

df_filtrado <- df_nuevo %>%
  # 2. Convertir de "24/3/2024" a Date
  mutate(ini_lact = dmy(ini_lact)) %>%
  # 3. Filtrar desde 1 de marzo de 2024 en adelante
  filter(ini_lact >= as.Date("2024-03-01"))
# Ver cuántos registros quedaron
nrow(df_filtrado)
## [1] 279

Renombramos las variables de interés para una mejor lectura.

abril_mayo <- abril_mayo %>%
  rename(
    grasa_may = Gr.UC,
    grasa_abr = Gr.CA,
    prot_may  = Pr.UC,
    prot_abr  = Pr.CA, 
    rp=RP
  )

Filtramos variables de interés, coincidentes con rp y renombramos para una mejor lectura

#Datos obtenidos del dataset abril_mayo
df_filtrado <- df_filtrado %>%
  left_join(
    abril_mayo %>% select(rp, grasa_may, grasa_abr, prot_may, prot_abr),
    by = "rp"
  )


#renombramos las variables de interés para una mejor lectura.
junio_julio <- junio_julio %>%
  rename(
    grasa_jun = Gr.CA,
    prot_jun  = Pr.CA, 
    grasa_jul = Gr.UC,
    prot_jul  = Pr.UC,
    rp=RP
  )

#Datos obtenidos del dataset junio_julio
df_filtrado <- df_filtrado %>%
  left_join(
    junio_julio %>% select(rp, grasa_jun, prot_jun, grasa_jul, prot_jul),
    by = "rp"
  )

Tratamiento de las variables de interés

#Convertir la fecha de inicio de lactancia a marzo, mayo o abril según corresponda para dejar valores correspondientes a mes 1 y mes 2 de lactancia.
df_filtrado <- df_filtrado %>%
  mutate(
    mes_lact = case_when(
      format(ini_lact, "%m") == "03" ~ "marzo",
      format(ini_lact, "%m") == "04" ~ "abril",
      format(ini_lact, "%m") == "05" ~ "mayo",
      TRUE                           ~ NA_character_
    )
  )


#Creamos las variables grasa_1 y grasa_2 correspondiente al valor de grasa obtenido en los dos meses posteriores al inicio de lactancia de la vaca
df_filtrado <- df_filtrado %>%
  mutate(
    grasa_1 = case_when(
      mes_lact == "marzo" ~ grasa_abr,
      mes_lact == "abril" ~ grasa_may,
      mes_lact == "mayo"  ~ grasa_jun,
      TRUE                ~ NA_real_
    ),
    grasa_2 = case_when(
      mes_lact == "marzo" ~ grasa_may,
      mes_lact == "abril" ~ grasa_jun,
      mes_lact == "mayo"  ~ grasa_jul,
      TRUE                ~ NA_real_
    )
  )

#Creamos las variables prot_1 y prot_2 correspondiente al valor de proteinas obtenidos en los dos meses posteriores al inicio de lactancia de la vaca
df_filtrado <- df_filtrado %>%
  mutate(
    prot_1 = case_when(
      mes_lact == "marzo" ~ prot_abr,
      mes_lact == "abril" ~ prot_may,
      mes_lact == "mayo"  ~ prot_jun,
      TRUE ~ NA_real_
    ),
    prot_2 = case_when(
      mes_lact == "marzo" ~ prot_may,
      mes_lact == "abril" ~ prot_jun,
      mes_lact == "mayo"  ~ prot_jul,
      TRUE ~ NA_real_
    )
  )


#Creamos las variables lts_1 y lts_2 correspondiente al valor de litros diario obtenidos en los dos meses posteriores al inicio de lactancia de la vaca

df_filtrado <- df_filtrado %>%
  mutate(
    rcs_1 = case_when(
      mes_lact == "marzo" ~ rcs_abr_24,
      mes_lact == "abril" ~ rcs_may_24,
      mes_lact == "mayo"  ~ rcs_jun_24,
      TRUE ~ NA_real_
    ),
    rcs_2 = case_when(
      mes_lact == "marzo" ~ rcs_may_24,
      mes_lact == "abril" ~ rcs_jun_24,
      mes_lact == "mayo"  ~ rcs_jul_24,
      TRUE ~ NA_real_
    )
  )



#Creamos las variables rcs_1 y rcs_2 correspondiente al valor de rcs obtenidos en los dos meses posteriores al inicio de lactancia de la vaca
df_filtrado <- df_filtrado %>%
  mutate(
    lts_1 = case_when(
      mes_lact == "marzo" ~ abr_24,
      mes_lact == "abril" ~ may_24,
      mes_lact == "mayo"  ~ jun_24,
      TRUE ~ NA_real_
    ),
    lts_2 = case_when(
      mes_lact == "marzo" ~ may_24,
      mes_lact == "abril" ~ jun_24,
      mes_lact == "mayo"  ~ jul_24,
      TRUE ~ NA_real_
    )
  )

# Unir con abril_mayo usando el nombre exacto de las columnas Vitales
df_filtrado <- df_filtrado %>%
  left_join(
    abril_mayo %>% 
      select(
        rp, 
        `Dias Vital`, 
        `Le.Vital`, 
        `Gr.Vital`, 
        `Pr.Vital`
      ),
    by = "rp"
  )
#Renombrar columnas 
junio_julio <- junio_julio %>%
  rename(
    diastotal_vital = `Dias Vital`,
    leche_vital     = Le.Vital,
    grasa_vital     = Gr.Vital,
    prot_vital      = Pr.Vital
  )


df_filtrado <- df_filtrado %>%
  left_join(
    junio_julio %>%
      select(
        rp,
        diastotal_vital,
        leche_vital,
        grasa_vital,
        prot_vital
      ),
    by = "rp"
  )

Creación de variables que expresan variación porcentual

#Agregamos variables que representan la variacion porcentual de un mes respecto de otro con el objetivo de que estas variables contengan la información para saber si hubo caída o subida de valores de un mes respecto de otro

df_filtrado <- df_filtrado %>%
  mutate(
    var_lts   = lts_2   - lts_1,
    var_grasa = grasa_2 - grasa_1,
    var_prot  = prot_2  - prot_1,
    var_rcs   = rcs_2   - rcs_1
  )



#Agregamos la cant de lts, grasa, prot y dias que acumula desde su inicio de lactancia hasta julio
df_filtrado <- df_filtrado %>%
  left_join(
    junio_julio %>%
      select(
        rp,
        leche_acum   = `Le.ACU`,
        grasa_acum = `Gr.ACU`,
        prot_acum  = `Pr.ACU`
      ),
    by = "rp"
  ) %>%
  mutate(
    dias_acum = as.integer(as.Date("2024-06-03") - ini_lact)
  )

Una vez obtenidas las variables predictoras, unimos todo en un mismo dataset

df_final <- df_filtrado %>%
  select(
    rp,
    grasa_1,
    grasa_2,
    prot_1,
    prot_2,
    lts_1,
    lts_2,
    rcs_1,
    rcs_2,
    lote_act,
    nro_lact, 
    le_305_acu, 
    var_lts, 
    var_grasa,
    var_prot,
    var_rcs,
    leche_acum,
    grasa_acum,
    prot_acum,
    dias_acum, 
    ini_lact,
    diastotal_vital,
    grasa_vital,
    prot_vital,
    leche_vital
  )

Diccionario del dataset en análisis

# Construir el diccionario
diccionario <- data.frame(
  variable_original = c(
    "rp", "grasa_1", "grasa_2", "prot_1", "prot_2",
    "lts_1", "lts_2", "rcs_1", "rcs_2", "lote_act",
    "nro_lact", "le_305_acu", "var_lts", "var_grasa", "var_prot", "var_rcs",
    "leche_acum", "grasa_acum", "prot_acum", "dias_acum",
    "diastotal_vital", "grasa_vital", "prot_vital", "leche_vital"
  ),
  nombre_variable = c(
    "identificador de la vaca",
    "grasa mes 1",
    "grasa mes 2",
    "proteína mes 1",
    "proteína mes 2",
    "litros día control mes 1",
    "litros día control mes 2",
    "RCS mes 1",
    "RCS mes 2",
    "lote actual",
    "número de lactancias",
    "leche acumulada 305 d",
    "∆ litros (mes 2-mes 1)",
    "∆ grasa (p.p.)",
    "∆ proteína (p.p.)",
    "∆ RCS",
    "leche acumulada",
    "grasa acumulada",
    "proteína acumulada",
    "días acumulados",
    "días totales vital",
    "grasa vital",
    "proteína vital",
    "leche vital"
  ),
  tipo = c(
    "categórico (chr)",
    rep("numérico (dbl)", 4),
    rep("numérico (dbl)", 2),
    rep("numérico (dbl)", 2),
    "categórico (chr)",
    "entero (int)",
    "numérico (dbl)",
    rep("numérico (dbl)", 4),
    rep("numérico (dbl)", 3),
    "entero (int)",
    "entero (int)",
    rep("numérico (dbl)", 3)
  ),
  descripcion = c(
    "número de identificación individual.",
    "porcentaje de grasa del 1.er mes (% p/v).",
    "porcentaje de grasa del 2.º mes (% p/v).",
    "porcentaje de proteína del 1.er mes (% p/v).",
    "porcentaje de proteína del 2.º mes (% p/v).",
    "litros producidos en el día de control del 1.er mes (L).",
    "litros producidos en el día de control del 2.º mes (L).",
    "recuento de células somáticas 1.er mes (×10³ cél/mL).",
    "recuento de células somáticas 2.º mes (×10³ cél/mL).",
    "grupo/lote productivo actual.",
    "número de lactancias (paridad).",
    "litros acumulados en los primeros 305 días de la lactancia actual (L).",
    "variación de litros entre mes 2 y mes 1 (L).",
    "variación de grasa entre mes 2 y mes 1 (puntos % ).",
    "variación de proteína entre mes 2 y mes 1 (puntos % ).",
    "variación de RCS entre mes 2 y mes 1 (×10³ cél/mL).",
    "litros acumulados en la lactancia actual (L).",
    "kg de grasa acumulados en la lactancia actual (kg).",
    "kg de proteína acumulados en la lactancia actual (kg).",
    "días transcurridos en la lactancia actual (días).",
    "días totales en producción desde el primer parto (días).",
    "kg de grasa producidos en toda la vida (kg).",
    "kg de proteína producidos en toda la vida (kg).",
    "litros de leche producidos en toda la vida (L)."
  ),
  stringsAsFactors = FALSE
)
# Mostrar el diccionario
kable(diccionario, "html", escape = FALSE, align = "l") %>% 
  kable_styling(full_width = FALSE,
                bootstrap_options = c("striped", "hover", "condensed"))
variable_original nombre_variable tipo descripcion
rp identificador de la vaca categórico (chr) número de identificación individual.
grasa_1 grasa mes 1 numérico (dbl) porcentaje de grasa del 1.er mes (% p/v).
grasa_2 grasa mes 2 numérico (dbl) porcentaje de grasa del 2.º mes (% p/v).
prot_1 proteína mes 1 numérico (dbl) porcentaje de proteína del 1.er mes (% p/v).
prot_2 proteína mes 2 numérico (dbl) porcentaje de proteína del 2.º mes (% p/v).
lts_1 litros día control mes 1 numérico (dbl) litros producidos en el día de control del 1.er mes (L).
lts_2 litros día control mes 2 numérico (dbl) litros producidos en el día de control del 2.º mes (L).
rcs_1 RCS mes 1 numérico (dbl) recuento de células somáticas 1.er mes (×10³ cél/mL).
rcs_2 RCS mes 2 numérico (dbl) recuento de células somáticas 2.º mes (×10³ cél/mL).
lote_act lote actual categórico (chr) grupo/lote productivo actual.
nro_lact número de lactancias entero (int) número de lactancias (paridad).
le_305_acu leche acumulada 305 d numérico (dbl) litros acumulados en los primeros 305 días de la lactancia actual (L).
var_lts ∆ litros (mes 2-mes 1) numérico (dbl) variación de litros entre mes 2 y mes 1 (L).
var_grasa ∆ grasa (p.p.) numérico (dbl) variación de grasa entre mes 2 y mes 1 (puntos % ).
var_prot ∆ proteína (p.p.) numérico (dbl) variación de proteína entre mes 2 y mes 1 (puntos % ).
var_rcs ∆ RCS numérico (dbl) variación de RCS entre mes 2 y mes 1 (×10³ cél/mL).
leche_acum leche acumulada numérico (dbl) litros acumulados en la lactancia actual (L).
grasa_acum grasa acumulada numérico (dbl) kg de grasa acumulados en la lactancia actual (kg).
prot_acum proteína acumulada numérico (dbl) kg de proteína acumulados en la lactancia actual (kg).
dias_acum días acumulados entero (int) días transcurridos en la lactancia actual (días).
diastotal_vital días totales vital entero (int) días totales en producción desde el primer parto (días).
grasa_vital grasa vital numérico (dbl) kg de grasa producidos en toda la vida (kg).
prot_vital proteína vital numérico (dbl) kg de proteína producidos en toda la vida (kg).
leche_vital leche vital numérico (dbl) litros de leche producidos en toda la vida (L).

Consideraciones sobre la variable leche vital y el número de parición

Se tuvo en cuenta que la variable leche vital representa la producción total de una vaca desde su primera lactancia. En los casos de primera parición, leche vital y leche acumulada fueron equivalentes, ya que no existía historial productivo anterior. En cambio, para vacas multíparas, leche vital incluyó la producción de ciclos previos.

Esta diferencia introdujo una variación estructural entre vacas primíparas y multíparas, lo cual podría haber afectado el comportamiento del modelo si no se contemplaba adecuadamente. Por tal motivo, se incorporó la variable número de pariciones como atributo en el entrenamiento, con el objetivo de que el modelo pudiera distinguir entre ambos perfiles y evitar sesgos en la interpretación de variables acumuladas.

Planteo del problema

Preguntas disparadoras

1- ¿Con solo 2 meses de datos podemos anticipar si la vaca superará 13 000 L? Flujo de predicción: datos Día 0-60 → Random Forest → proyección Día 305

Se evalúa si los registros productivos de los primeros 60 días de lactancia permiten anticipar, con suficiente precisión, si una vaca superará los 13 000 L al día 305.

2- ¿Para qué predecir algo que ya conocemos? Utilizamos datos históricos para entrenar al modelo, esto permite aprender patrones y que el modelo aprenda.

Idea central

El algoritmo Random Forest utiliza únicamente esa “foto temprana” de dos meses para proyectar el rendimiento al final de la lactancia.

Entrenamiento del modelo

  • Conjunto de datos: historiales completos de vacas del propio tambo que ya finalizaron su lactancia.
  • Objetivo: aprender los patrones de arranque característicos de vacas de alta versus baja producción.

Aplicaciones prácticas

  1. Próxima lactancia de vacas existentes: permite ajustar dieta, manejo o plantear el descarte con antelación.
  2. Ingreso de vacas nuevas: facilita decidir la continuidad o la inversión sin esperar un ciclo completo.

Ventajas de utilizar Random Forest

  • Robusto frente a múltiples variables y relaciones no lineales.
  • Proporciona un ranking de importancia de los predictores, útil para priorizar acciones de manejo.

Analogía explicativa

Así como un docente proyecta el desempeño anual a partir del primer trimestre, el modelo pronostica toda la lactancia basándose en las primeras ocho semanas.

Conclusión

No se analizan los datos históricos por curiosidad, sino para informar decisiones futuras y maximizar el valor de la información que el tambo genera de forma rutinaria.

Análisis breve del nuevo dataset

# 1. Contar NAs por columna
na_counts <- colSums(is.na(df_final))
print(na_counts)
##              rp         grasa_1         grasa_2          prot_1          prot_2 
##               0               1               0               1               0 
##           lts_1           lts_2           rcs_1           rcs_2        lote_act 
##               1               0               1               0               0 
##        nro_lact      le_305_acu         var_lts       var_grasa        var_prot 
##               0               0               1               1               1 
##         var_rcs      leche_acum      grasa_acum       prot_acum       dias_acum 
##               1               0               0               0               0 
##        ini_lact diastotal_vital     grasa_vital      prot_vital     leche_vital 
##               0               0               0               0               0
#Los datos faltantes se deben a un error del sistema, se imputan con la media.
df_final <- df_final %>%
  mutate(across(
    c(grasa_1, prot_1, lts_1, rcs_1,
      var_lts, var_grasa, var_prot, var_rcs),
    ~ ifelse(is.na(.x), mean(.x, na.rm = TRUE), .x)
  ))
# Verificar que ya no queden NAs
colSums(is.na(df_final))
##              rp         grasa_1         grasa_2          prot_1          prot_2 
##               0               0               0               0               0 
##           lts_1           lts_2           rcs_1           rcs_2        lote_act 
##               0               0               0               0               0 
##        nro_lact      le_305_acu         var_lts       var_grasa        var_prot 
##               0               0               0               0               0 
##         var_rcs      leche_acum      grasa_acum       prot_acum       dias_acum 
##               0               0               0               0               0 
##        ini_lact diastotal_vital     grasa_vital      prot_vital     leche_vital 
##               0               0               0               0               0
#Estadísticas básicas de todo el data.frame
summary(df_final)
##        rp           grasa_1         grasa_2          prot_1     
##  Min.   :14012   Min.   :1.550   Min.   :1.130   Min.   :2.640  
##  1st Qu.:16570   1st Qu.:2.940   1st Qu.:2.210   1st Qu.:3.330  
##  Median :17340   Median :3.800   Median :2.760   Median :3.550  
##  Mean   :17374   Mean   :3.871   Mean   :2.742   Mean   :3.750  
##  3rd Qu.:18549   3rd Qu.:4.545   3rd Qu.:3.215   3rd Qu.:3.935  
##  Max.   :19507   Max.   :9.570   Max.   :5.960   Max.   :7.160  
##      prot_2          lts_1           lts_2           rcs_1       
##  Min.   :2.550   Min.   :15.40   Min.   :32.40   Min.   :   0.0  
##  1st Qu.:3.105   1st Qu.:36.70   1st Qu.:43.90   1st Qu.:  12.0  
##  Median :3.300   Median :43.40   Median :50.60   Median :  37.0  
##  Mean   :3.295   Mean   :43.12   Mean   :51.36   Mean   : 168.3  
##  3rd Qu.:3.450   3rd Qu.:49.60   3rd Qu.:58.50   3rd Qu.: 137.5  
##  Max.   :4.240   Max.   :67.00   Max.   :76.00   Max.   :2916.0  
##      rcs_2           lote_act        nro_lact       le_305_acu   
##  Min.   :   1.0   Min.   :3.000   Min.   :1.000   Min.   : 7755  
##  1st Qu.:   5.0   1st Qu.:4.000   1st Qu.:1.000   1st Qu.:12773  
##  Median :  11.0   Median :5.000   Median :2.000   Median :14273  
##  Mean   : 104.1   Mean   :5.219   Mean   :1.781   Mean   :14130  
##  3rd Qu.:  22.5   3rd Qu.:6.000   3rd Qu.:2.000   3rd Qu.:15830  
##  Max.   :5107.0   Max.   :7.000   Max.   :5.000   Max.   :18698  
##     var_lts          var_grasa         var_prot          var_rcs        
##  Min.   :-12.400   Min.   :-5.990   Min.   :-3.7200   Min.   :-2791.00  
##  1st Qu.:  3.000   1st Qu.:-1.825   1st Qu.:-0.5950   1st Qu.:  -90.00  
##  Median :  7.000   Median :-0.990   Median :-0.2400   Median :  -16.00  
##  Mean   :  8.266   Mean   :-1.130   Mean   :-0.4559   Mean   :  -63.85  
##  3rd Qu.: 12.500   3rd Qu.:-0.300   3rd Qu.:-0.0350   3rd Qu.:   -1.00  
##  Max.   : 42.400   Max.   : 2.140   Max.   : 0.5800   Max.   : 5091.00  
##    leche_acum     grasa_acum      prot_acum        dias_acum    
##  Min.   :1961   Min.   : 41.0   Min.   : 67.32   Min.   :26.00  
##  1st Qu.:4027   1st Qu.:117.3   1st Qu.:136.97   1st Qu.:53.50  
##  Median :4920   Median :144.2   Median :169.73   Median :75.00  
##  Mean   :4856   Mean   :142.1   Mean   :164.89   Mean   :67.51  
##  3rd Qu.:5773   3rd Qu.:169.0   3rd Qu.:195.54   3rd Qu.:84.00  
##  Max.   :8036   Max.   :247.6   Max.   :267.37   Max.   :94.00  
##     ini_lact          diastotal_vital   grasa_vital       prot_vital     
##  Min.   :2024-03-01   Min.   :  59.0   Min.   :  44.2   Min.   :  67.32  
##  1st Qu.:2024-03-11   1st Qu.: 113.0   1st Qu.: 146.2   1st Qu.: 164.56  
##  Median :2024-03-20   Median : 411.0   Median : 577.8   Median : 611.03  
##  Mean   :2024-03-27   Mean   : 374.6   Mean   : 572.0   Mean   : 555.17  
##  3rd Qu.:2024-04-10   3rd Qu.: 562.0   3rd Qu.: 873.9   3rd Qu.: 806.77  
##  Max.   :2024-05-08   Max.   :1317.0   Max.   :2158.0   Max.   :2060.57  
##   leche_vital   
##  Min.   : 1961  
##  1st Qu.: 4684  
##  Median :16936  
##  Mean   :16288  
##  3rd Qu.:22840  
##  Max.   :63588
#Definimos la variable OBJETIVO
df_final <- df_final %>%
  mutate(
    condicion = ifelse(le_305_acu > 13000, "alta", "baja")
  )
# Verifica distribución de la nueva variable
df_final %>%
  count(condicion)
## # A tibble: 2 × 2
##   condicion     n
##   <chr>     <int>
## 1 alta        201
## 2 baja         78

6. Exploración inicial de datos (EDA) en dataset final

Distribución de la variable objetivo “condición” que define alta o baja producción de las vacas.

# Paso 1: Distribución de la variable objetivo "condicion"
# a) Cálculo de frecuencias y proporciones
freq_target <- table(df_final$condicion)
prop_target <- prop.table(freq_target)

Proporción en cantidad y %

print(freq_target)               # Conteo de vacas por clase (baja / alta)
## 
## alta baja 
##  201   78
print(round(prop_target, 3))     # Proporciones redondeadas a 3 decimales
## 
## alta baja 
## 0.72 0.28

Gráfico ilustrativo de la proporción

# b) Gráfico de barras con ggplot2 para visualizar el balance de clases
if (!requireNamespace("ggplot2", quietly = TRUE)) {
  install.packages("ggplot2")
}

ggplot(df_final, aes(x = condicion)) +
  geom_bar(fill = "steelblue") +
  labs(
    title = "Distribución de la producción: baja vs alta",
    x = "Clase de producción",
    y = "Número de vacas"
  ) +
  theme_minimal()

En nuestro set de datos hay muchas más vacas de “alta” producción (72%) que de “baja” (28%). A primera vista eso podría parecer un sesgo, pero como nuestro objetivo es precisamente detectar bien las vacas de alta producción, ese desequilibrio juega a nuestro favor: un Random Forest entrenado con esta misma proporción aprende de manera natural a reconocer esa clase mayoritaria con un recall alto.

Intentar equilibrar el dataset trae más problemas que soluciones:

  • SMOTE u otros oversampling crea vacas “falsas” que no reflejan variabilidad real.

  • Undersampling elimina datos verdaderos de alta producción, quitando riqueza al entrenamiento.

Por eso, en lugar de generar o quitar casos, podemos confiar en la distribución natural y, si más adelante quisiéramos mejorar aún más la sensibilidad, ajustaríamos el umbral de clasificación o aplicaríamos pesos al evaluar el modelo.

Análisis univariado de las variables en estudio

# Función para plotear un grupo de variables
plot_counts <- function(vars, title) {
  df_final %>% 
    select(all_of(vars)) %>% 
    pivot_longer(cols = everything(), names_to = "variable", values_to = "value") %>% 
    ggplot(aes(x = value)) +
    geom_histogram(bins = 20, fill = "steelblue", alpha = 0.7) +
    facet_wrap(~ variable, scales = "free", ncol = min(3, length(vars))) +
    labs(title = title, x = NULL, y = "Cantidad") +
    theme_minimal()
}

Grasa

plot_counts(
  c("grasa_1", "grasa_2", "var_grasa", "grasa_acum", "grasa_vital"),
  "Grasa: mes 1, mes 2, variación, acumulado y vital"
)

– En el primer mes (grasa_1) la mayoría se mueve entre 3–6 % y en el segundo (grasa_2) entre 2–4 %, con ligera asimetría a la derecha.

– La variación (var_grasa) es mayoritariamente negativa (caída de 1–3 %), lo que indica un descenso de la grasa promedio del mes 1 al 2.

– El acumulado (grasa_acum) y el vital (grasa_vital) muestran distribuciones en forma de campana con cola a la derecha, señalando unos pocos animales con acumulados muy altos.

Proteína

plot_counts(
  c("prot_1", "prot_2", "var_prot", "prot_acum", "prot_vital"),
  "Proteína: mes 1, mes 2, variación, acumulado y vital"
)

– En el primer mes (prot_1), la mayoría de vacas se ubica entre 3.3–3.9 %, con mínimos en 2.64 % y algunos picos aislados hasta 7.16 %, lo que genera una cola derecha leve.

– En el segundo mes (prot_2), la concentración es más estrecha, entre 3.1–3.45 %, con escasas observaciones por encima de 4 % y ligera asimetría hacia valores mayores.

– La variación (var_prot) es negativa para la mayoría (caídas de 0.24–0.6 %), con unos pocos outliers de descenso fuerte y muy pocos incrementos positivos (hasta +0.58 %).

– El acumulado (prot_acum) forma una campana con cola derecha, centrado en 137–196 unidades, señalando varios animales con producción total de proteína destacada.

– El vital (prot_vital) muestra fuerte sesgo derecho: gran parte de los valores entre 165–807, pero varios outliers por encima de 1 500.

Cantidad de litros de leche

plot_counts(
  c("lts_1", "lts_2", "var_lts", "leche_acum", "leche_vital"),
  "Leche: mes 1, mes 2, variación, acumulado y vital"
)

– lts_1 y lts_2 concentran la mayoría entre 35–60 L/día, con forma unimodal y ligera asimetría a la derecha.

– La variación (var_lts) es mayoritariamente positiva (incrementos de ~5–20 L), reflejando que muchas vacas aumentan su producción en el segundo mes.

– El acumulado en dos meses (leche_acum) se agrupa en torno a 3 000–7 000 L, con unos pocos ejemplares por encima de ese rango.

– La producción vital (leche_vital) muestra fuerte sesgo derecho, con la mayoría entre 20 000–40 000 L y varios outliers hasta ~60 000 L.

Cantidad de días

plot_counts(
  c("dias_acum", "diastotal_vital"),
  "Días: acumulados y totales de vida productiva"
)

– dias_acum (primeros dos meses) se concentra entre 75–90 días, con algunos valores por debajo de 50 días.

– diastotal_vital (duración total) muestra un grupo fuerte en 100–200 días (periodos cortos de lactancia) y otro amplio entre 300–700 días, con unos pocos animales superando los 1 000 días, lo que puede indicar lactancias muy largas.

Lote y número de lactancia

plot_counts(
  c("lote_act", "nro_lact"),
  "Contexto: lote actual y número de lactancias previas"
)

– lote_act agrupa la mayoría de vacas en los lotes 4–7, con un pico en el lote 4 (~105 vacas) y muy pocos en el lote 3 (~10 vacas). Esto aporta información sobre el manejo por grupos homogéneos.

– nro_lact revela que la mayor parte de los animales está en su 1.ª o 2.ª lactancia (~120 y ~110 vacas), cayendo drásticamente a partir de la tercera lactancia. Hay muy pocos casos con ≥4 lactancias, por lo que el modelo contará con escasa señal para vacas de larga trayectoria.

Recuento de Células Somáticas (RCS)

plot_counts(
  c("rcs_1", "rcs_2"),
  "Recuento de células somáticas: mes 1 y mes 2"
)

– En ambos meses (rcs_1, rcs_2), la mayoría de vacas registra recuentos bajos (<100), y se observa una fuerte asimetría derecha con unos pocos valores extremos (hasta >3 000 en mes 1 y >4 000 en mes 2).

– Estos outliers pueden reflejar problemas de salud mamaria.

7. Aplicación del modelo Random Forest

1. Conversión de la variable “condicion” a factor

# Convertir a factor
df_final <- df_final %>%
  mutate(condicion = factor(condicion, levels = c("baja", "alta")))
# Verificar estructura y frecuencia
str(df_final$condicion)
##  Factor w/ 2 levels "baja","alta": 2 2 2 2 2 2 1 2 2 1 ...
table(df_final$condicion)
## 
## baja alta 
##   78  201

2. Partición entrenamiento/prueba (70/30 estratificado)

set.seed(123)
train_idx  <- createDataPartition(df_final$condicion, p = 0.7, list = FALSE)
train_data <- df_final[train_idx, ]
test_data  <- df_final[-train_idx, ]

Ver tamaño de cada conjunto

Entrenamiento:

print(table(train_data$condicion))
## 
## baja alta 
##   55  141

Prueba:

print(table(test_data$condicion))
## 
## baja alta 
##   23   60

3. Entrenar el Random Forest

set.seed(123)
modelo_rf <- randomForest(
  condicion ~ grasa_1 + grasa_2 + prot_1 + prot_2 +
    lts_1   + lts_2   + rcs_1   + rcs_2 +
    var_lts + var_grasa + var_prot + var_rcs +
    leche_acum + grasa_acum + prot_acum + dias_acum +
    lote_act + nro_lact+ diastotal_vital + grasa_vital +  prot_vital + leche_vital,
  data       = train_data,
  ntree      = 200,
  mtry       = 4,          
  importance = TRUE
)

Elección de parámetros ntree=200, es un número razonable: evita sobreentrenamiento, pero da buena estabilidad y precisión.** mtry=4, se suele usar sqrt(p) para clasificación, donde p es el número total de predictores, tenemos 22 variables → sqrt(22) ≈ 4.69 → se usó 4.

# Mostrar resumen del modelo
print(modelo_rf)
## 
## Call:
##  randomForest(formula = condicion ~ grasa_1 + grasa_2 + prot_1 +      prot_2 + lts_1 + lts_2 + rcs_1 + rcs_2 + var_lts + var_grasa +      var_prot + var_rcs + leche_acum + grasa_acum + prot_acum +      dias_acum + lote_act + nro_lact + diastotal_vital + grasa_vital +      prot_vital + leche_vital, data = train_data, ntree = 200,      mtry = 4, importance = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 200
## No. of variables tried at each split: 4
## 
##         OOB estimate of  error rate: 11.73%
## Confusion matrix:
##      baja alta class.error
## baja   42   13  0.23636364
## alta   10  131  0.07092199

Al evaluar el desempeño interno (OOB error ~12%), veo que el modelo acierta en torno al 88% de los casos y, lo más importante, detecta casi todas las vacas de alta producción con apenas un ~7% de error. Dado que nuestro objetivo es identificar aquellas de alto rendimiento, esto indica que rara vez se nos escapará una vaca buena. No obstante, existe un pequeño porcentaje de vacas de baja producción que puede recibir la etiqueta “alta”; si resulta crítico minimizar falsas expectativas, podríamos ajustar el umbral de decisión o incorporar un criterio adicional antes de clasificar definitivamente como “alta”.

4. Evaluación en el conjunto de prueba

# Predecir en test
pred_test <- predict(modelo_rf, newdata = test_data)
# Matriz de confusión y métricas
conf_mat <- confusionMatrix(pred_test, test_data$condicion, positive="alta")
print(conf_mat)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction baja alta
##       baja   18    4
##       alta    5   56
##                                           
##                Accuracy : 0.8916          
##                  95% CI : (0.8041, 0.9492)
##     No Information Rate : 0.7229          
##     P-Value [Acc > NIR] : 0.0001681       
##                                           
##                   Kappa : 0.7257          
##                                           
##  Mcnemar's Test P-Value : 1.0000000       
##                                           
##             Sensitivity : 0.9333          
##             Specificity : 0.7826          
##          Pos Pred Value : 0.9180          
##          Neg Pred Value : 0.8182          
##              Prevalence : 0.7229          
##          Detection Rate : 0.6747          
##    Detection Prevalence : 0.7349          
##       Balanced Accuracy : 0.8580          
##                                           
##        'Positive' Class : alta            
## 

Matriz de confusión

  • Al ver que el modelo predice correctamente 56 vacas “alta” y 18 vacas “baja”, y se equivoca en 4 falsos negativos y 5 falsos positivos, entiendo dónde concentran mis aciertos y errores: acierto más al etiquetar “alta” y cometo pocos fallos al descartar “baja”.

Accuracy (89.16 %)

  • Casi 9 de cada 10 predicciones son correctas. Esto me da confianza en el rendimiento general, aunque sé que debo mirar más métricas para entender bien los desequilibrios.

IC 95 % (0.8041–0.9492)

  • Estoy 95 % seguro de que la verdadera accuracy, si repitiera el experimento, quedaría entre 80,4 % y 94,9 %. Ese margen me muestra que mi estimación es estable.

No Information Rate (72.29 %)

  • Si siempre dijera “alta” (la clase más frecuente), acertaría el 72,3 %. Mi modelo supera con creces esa línea de base, lo que indica un verdadero valor predictivo.

P-Value [Acc > NIR] (0.00017)

  • La probabilidad de obtener mi accuracy (89,16 %) solo por azar, asumiendo que la mejor estrategia fuera predecir siempre “alta”, es prácticamente cero. Esto confirma que mi modelo supera significativamente la regla más tonta.

Kappa (0.7257)

  • Un Kappa de 0,73 revela un acuerdo fuerte entre predicción y realidad, más allá de lo explicable por la frecuencia de cada clase. Me asegura que no estoy “adivinando” con base en la prevalencia.

Mcnemar’s Test P-Value (1.00)

  • Con p=1, no hay sesgo sistemático en los errores: no tiendo a confundir siempre “alta” con “baja” o viceversa. Mis falsos positivos y negativos están equilibrados.

Sensibilidad (93.33 %)

  • De todas las vacas que realmente son de alta producción, el modelo identifica el 93.3 %. Esto es clave: casi no dejo escapar ejemplares valiosos.

Especificidad (78.26 %)

  • De todas las vacas de baja producción, el 78.3 % las detecto correctamente. Aunque es menor que la sensibilidad, sigue siendo un nivel sólido de aciertos para la clase minoritaria.

Precisión / PPV (91.80 %)

  • Cuando el modelo dice “alta”, lo hace bien el 91.8 % de las veces. Puedo confiar en que pocas de mis señales positivas son falsas alarmas.

NPV (81.82 %)

  • Cuando predigo “baja”, acierto el 81.8 % de las veces. Esto me ayuda a no descartar vacas que podrían ser buenas candidatas.

Prevalencia (72.29 %)

  • Refleja que el 72.3 % de las vacas en mi muestra son realmente “alta”. Tomar en cuenta esta proporción me evita sobreestimar el rendimiento del modelo.

Detection Rate (67.47 %)

  • El 67.5 % de todas las vacas (prevalencia × sensibilidad) son correctamente detectadas como “alta”. Es una métrica combinada que muestra el alcance real de la clase positiva.

Detection Prevalence (73.49 %)

  • El 73.5 % de las predicciones son “alta” (incluye falsos positivos). Me indica cuán generoso es el modelo al asignar la etiqueta positiva.

Balanced Accuracy (85.80 %)

  • Promedia sensibilidad y especificidad, dándome un 85.8 %. Este valor ponderado me asegura que desempeño bien con ambas clases, incluso frente al desbalance.

‘Positive’ Class: alta

  • Todas estas métricas de sensibilidad y precisión se refieren a la clase “alta”, que es mi principal interés para maximizar la identificación de vacas de alto rendimiento.

Rendimiento global y confianza estadística

Al observar la exactitud en el conjunto de prueba, noto que el modelo acierta alrededor de 0.8916 (≈89.2%) de las predicciones, y el intervalo de confianza al 95% va de 0.804 a 0.949, lo que indica que, aun considerando variabilidad muestral, la precisión se mantiene en niveles elevados. Comparado con la tasa de acierto si siempre predijéramos la clase mayoritaria (“alta”), que sería de ≈72.3%, el p-valor asociado (≈0.00017) confirma que el rendimiento real supera significativamente esa línea base sin ser producto del azar. El coeficiente Kappa (~0.726) refuerza que el acuerdo entre predicción y realidad va más allá de lo que correspondería por la distribución de clases. En conjunto, estos indicadores me hacen confiar en que el modelo captura señales relevantes de las variables y no se limita a reflejar la prevalencia de la clase más frecuente.

Desglose por clases y enfoque en alta producción

El modelo identifica correctamente casi 8 de cada 10 vacas “alta” (sensibilidad ≈78.3%) y, de las predicciones “alta”, el 81.8% son efectivamente buenas, lo cual indica que la mayoría de etiquetas positivas realmente corresponden a animales valiosos. Además, las vacas clasificadas como “baja” son correctas en un 91.8%, por lo que la confiabilidad del modelo en ambas clases es elevada y se pueden tomar decisiones bastante seguras a partir de sus resultados.

Equilibrio de sensibilidad y especificidad, y distribución de clases

Aunque “alta” representa la mayoría del total de vacas (≈72%), el modelo logra equilibrar bien ambas clases: detecta con buena sensibilidad a las vacas “alta” y mantiene un excelente desempeño para identificar correctamente las “baja”, alcanzando una balanced accuracy de ≈0.858. El test de McNemar no evidencia sesgo sistemático en los errores, lo cual sugiere que no hay una inclinación significativa a favorecer una clase por sobre otra en los desaciertos.

Implicaciones prácticas y recomendaciones

Para el tambo, el modelo Random Forest resulta muy útil para detectar con alta confianza vacas de alto rendimiento, lo que permite orientar recursos, manejos y decisiones reproductivas hacia animales de mayor valor. Si bien hay un pequeño margen de error en ambas clases (falsos positivos y negativos), el nivel de acierto permite confiar en la herramienta como una primera capa de decisión. Para decisiones muy críticas (descarte o selección genética), se recomienda complementarlo con monitoreo adicional o umbrales conservadores. En términos generales, el modelo ofrece una base sólida y confiable para identificar ejemplares con potencial lechero destacado.

5. Importancia de variables

# Extraer importancia numérica
imp <- importance(modelo_rf)
print(imp)
##                         baja        alta MeanDecreaseAccuracy MeanDecreaseGini
## grasa_1         -1.701450084  1.31956431            0.1339347        1.5137142
## grasa_2          0.694797704  1.44515802            1.5556654        1.7796242
## prot_1           0.440752039  1.89851183            1.9492854        1.4494098
## prot_2           4.293781223  1.10925027            3.9370951        3.4053819
## lts_1            7.674018470  6.32079986            9.9462125        6.8847189
## lts_2           10.298311330 10.73217473           13.5377294       13.3693194
## rcs_1           -0.352064257  2.09699608            1.4547093        1.5092482
## rcs_2           -0.420596081  0.98237033            0.4613079        1.2306780
## var_lts          0.279892211  2.79340280            2.5101306        1.7812422
## var_grasa        0.934926085  1.60572702            2.0578310        1.7509788
## var_prot         1.519523936  2.37451277            2.7119514        1.6541942
## var_rcs          0.078335736  2.85770413            2.3751337        1.5375804
## leche_acum       8.143369620  7.69895756           10.0640921        8.6234579
## grasa_acum       0.006969416  2.25999782            1.8429002        1.9581162
## prot_acum        3.208693637  4.47755245            5.8749013        4.4194652
## dias_acum       -0.306245131  2.82674362            2.1512193        1.7124133
## lote_act         1.426198411 -0.03086936            0.7471052        0.8841565
## nro_lact         0.295912239  2.56299529            2.6658905        1.0369271
## diastotal_vital  2.267982561  4.94214657            5.0913054        2.6748508
## grasa_vital     -0.483999084  4.73497458            4.2865588        3.1779634
## prot_vital       3.934244873  4.96531289            7.2022463        5.7416586
## leche_vital      8.761152044  7.85612156           10.7326647       10.8025539
# Gráfico de las Top 10 variables
varImpPlot(modelo_rf, n.var = 10, main = "Top 10 Variables Importantes")

Las variables más influyentes describen cuánta leche y cuántos sólidos produce la vaca muy temprano y qué tan constante ha sido a lo largo de su vida:

  • lts_2 (día 60): un volumen alto a los dos meses indica buena salud y un balance energético favorable; las vacas que ya producen mucho a esa altura suelen mantener rendimientos elevados el resto del ciclo.

  • lts_1 (día 30) y leche_acum (litros acumulados al corte): muestran la rapidez y la consistencia con que la vaca inicia la lactancia; un arranque firme anticipa mayores litros totales.

  • leche_vital: resume genética, manejo y sanidad de toda la vida productiva; quienes acumularon más litros en lactancias anteriores tienden a repetir ese comportamiento.

  • prot_acum / prot_vital: la proteína es un sólido estable y ligado a eficiencia alimentaria; valores altos refuerzan la señal de una vaca productiva y de buena calidad de leche.

  • diastotal_vital: contar con muchos días útiles de ordeñe refleja longevidad y ausencia de problemas reproductivos o mamarios; eso se traduce en más litros por ciclo.

  • grasa_vital: complementa al volumen; litros con sólidos elevados aportan mejor retorno económico.

En contraste, variables como RCS o el porcentaje de grasa del primer control muestran mayor variabilidad y menor poder explicativo en este conjunto de datos. En síntesis, los litros y sólidos medidos en los primeros 60 días, junto con el historial productivo, concentran la información clave para anticipar la producción al día 305.

Desde una perspectiva operativa, estos resultados sugieren que conviene centrar los esfuerzos de registro y validación de datos en las variables de volumen de producción y proteínas, optimizando así la eficiencia del modelo. Reducir o ajustar aquellas variables con escaso aporte puede ayudar a simplificar el modelo, evitar sobreajuste y facilitar su aplicación práctica en la gestión del rodeo.

8. Fundamentación del valor del modelo a través de la variabilidad mensual

Medir la varianza y el desvío estándar de la producción mensual nos permite cuantificar cuánta fluctuación tiene cada vaca mes a mes.

  • Si la variabilidad es muy baja, basta un simple corte univariable para predecir su rendimiento.

  • Si es alta, un modelo multivariable como Random Forest, que integra litros, grasa, proteína y RCS, añade valor real al capturar patrones complejos.

varianza <- df_modificado %>%
  # Filtrar solo los rp que aparecen en df_filtrado
  filter(rp %in% df_filtrado$rp) %>%
  # Seleccionar las columnas deseadas
  select(
    rp,
    abr_24,
    may_24,
    jun_24,
    jul_24,
    ago_24,
    sept_24,
    oct_24,
    nov_24,
    dic_24,
    ene_25,
    feb_25,
    mar_25
  )


varianza <- varianza %>%
  mutate(
    n_missing = rowSums(
      is.na(
        select(., abr_24, may_24, jun_24, jul_24, ago_24, sept_24,
               oct_24, nov_24, dic_24, ene_25, feb_25, mar_25)
      )
    )
  )


#cálculo de varianza de cada registro
varianza$var_meses <- apply(
  varianza[, c("abr_24","may_24","jun_24","jul_24","ago_24","sept_24",
               "oct_24","nov_24","dic_24","ene_25","feb_25","mar_25")],
  1,
  var,
  na.rm = TRUE
)
print(varianza$var_meses)
##   [1] 224.912727 212.460606 130.643636 102.682182 240.924545  56.806545
##   [7]  77.163636 216.066182 315.416000  73.038545  75.194182  85.989091
##  [13]  58.152727  92.128727 179.842182  50.556061 105.954909 276.931636
##  [19] 155.904444 140.450909 144.218182  78.859636  48.919556  85.762222
##  [25]  63.678545 144.044545  65.159697 112.406545 272.826909  38.888727
##  [31] 101.015556 108.360000 177.779636 114.298909 118.905455  99.239697
##  [37] 257.831515 172.834182  54.192424  71.130182  48.619636 322.810909
##  [43] 193.968444 167.737333 116.354909  40.658182 121.886061 188.835636
##  [49] 157.996061  80.999697 273.138909 227.708444 194.218909  64.049697
##  [55] 238.588889  66.130182  31.368727 126.753778 137.414242  16.517778
##  [61]  57.166970 379.315152  44.698909  64.198788 175.194909 151.211636
##  [67]  33.792727  86.250909  56.779091 186.890606  20.296000  53.216970
##  [73] 173.666364  65.710000 202.448727 118.009697 294.442909 150.971636
##  [79]  66.160727  63.536970  76.679697  75.510545  42.904727 173.137778
##  [85] 123.619636 343.410909  47.610606 105.078545 106.526222  50.944242
##  [91]  42.958545 135.440909  86.949091  91.202182 200.491636  49.251556
##  [97]  84.140606 109.695455  88.602909  58.621515  52.887273 216.110545
## [103] 156.488182  41.246970  71.853333  95.062545  35.904242 100.688000
## [109]  51.516061  88.439697  48.684889 154.528727 135.377879  75.002182
## [115]  70.019636  48.351515 171.240000  40.716061  15.676061 321.710667
## [121]  70.305455  80.671515 224.046061  51.735379 180.970667 298.803333
## [127]  48.134242  50.883333  23.418182  14.175455  27.518545  67.724889
## [133]  66.531818  19.734242  90.451818  37.695455 148.003636  48.658182
## [139]  91.993636  76.826667  12.308788  35.179091  80.131636  69.040909
## [145] 232.616970  73.715152  11.021515  64.576727  58.985152  24.355152
## [151] 184.134242  27.771636  74.002424  53.280444  42.784242  67.082424
## [157]  25.602424 129.856000 223.283333  85.102667  75.230000  92.636061
## [163]  30.064242  48.662727  43.370606  10.342545  29.602424  24.533333
## [169]  35.030545  24.720000  35.486970  55.217879  60.715152  57.317273
## [175] 109.867879  58.569697   9.900444  18.017879  51.672727  23.219636
## [181]  71.145152  31.175111 188.530909 142.017879  42.629091  25.202424
## [187]  67.814242  35.553778 113.464242  43.368727  29.866909  47.665455
## [193]   7.450909 119.782545  56.678545  68.461515  38.069091 122.366970
## [199]  13.653333  35.007273  37.075152  14.022727 137.544000  16.265152
## [205]  33.120909 132.997273  60.498909 140.771818  22.621818  18.831515
## [211]  10.294242  14.176000  40.499556  88.730606  10.073636  63.478788
## [217]  20.779091  22.570909  14.090909  15.564545  30.784727  89.950000
## [223]  33.380606  13.792424  15.566667  22.211818  36.737879  77.283333
## [229]  20.834909  64.578182  12.024242  49.474182  65.933333   9.447273
## [235]  74.519556  81.625455  48.522909 187.033636  18.946182  82.708788
## [241]  76.432424  10.250606  26.322182  40.794909  90.182545  42.774242
## [247] 142.152000   8.444889  42.042182  55.120727  62.795636 102.206667
## [253]  46.597273  36.931818  34.451556  47.240000  82.402909  97.334545
## [259]  12.231556  16.614242  10.308788  62.738182  30.642222  14.174333
## [265]   8.791556  69.317333  82.307636  48.021778  29.442424 156.417778
## [271] 110.997778  71.023222  16.149333  53.758222  16.110667  22.907111
## [277]  37.923636  10.551556  38.571111
# Calcula el promedio de la columna 'var_meses'
promedio_varianza <- mean(varianza$var_meses, na.rm = TRUE)
# Muestra el resultado
print(promedio_varianza)
## [1] 85.82655

cálculo del desvío estándar

Calculamos el desvío estándar porque, al medirse en litros, refleja de forma directa y comprensible la variabilidad mensual de cada vaca, permite fijar umbrales claros (p. ej. >10 L) y, así, decidir objetivamente cuándo un modelo multivariable aporta verdadero valor frente a una regla basada en un solo mes.

# Calcular desviación estándar de la producción mensual por vaca
varianza <- varianza %>%
  mutate(
    sd_meses = apply(
      select(., abr_24, may_24, jun_24, jul_24, ago_24, sept_24,
             oct_24, nov_24, dic_24, ene_25, feb_25, mar_25),
      1,
      sd,
      na.rm = TRUE
    )
  )

Calculamos el percentil 75 de la desviación estándar con estos objetivos:

1- Definir un umbral cuantitativo: obtener un valor robusto (no afectado por outliers) que nos permita separar las vacas con variabilidad “alta” de las de variabilidad “baja”.

2- Identificar el subgrupo que más necesita el modelo multivariable: las vacas por encima de Q3 (>10.63 L) presentan fluctuaciones fuertes mes a mes, de modo que no basta un único punto de corte y ahí el Random Forest añade mayor precisión.

3- Orientar la interpretación y el uso práctico: distinguir cuándo basta con una regla simple (sd baja) o cuándo conviene aplicar plenamente el modelo para reducir errores en casos de alta oscilación.

# Calcular el percentil 75 de la desviación estándar
q3_sd <- quantile(varianza$sd_meses, probs = 0.75, na.rm = TRUE)
print(q3_sd)
##      75% 
## 10.62708
# Etiquetar vacas con desviación estándar alta
varianza <- varianza %>%
  mutate(
    sd_alta = sd_meses > q3_sd
  )

# Calcular porcentaje de vacas con desviación alta
varianza %>%
  summarise(
    pct_sd_alta = mean(sd_alta, na.rm = TRUE) * 100
  )
## # A tibble: 1 × 1
##   pct_sd_alta
##         <dbl>
## 1        25.1

Este 25,1 % (unos 70 de 279 animales) corresponde a las vacas cuya producción fluctúa más de 10,6 L mes a mes. En la práctica, esto significa:

1- Si aplicáramos solo un criterio univariable (“si mes 1 > umbral → alta”), en este 25 % de vacas tendríamos una tasa de errores sensiblemente mayor:

  • Falsos negativos (descartar ejemplares de alto potencial) implican pérdidas de producción durante toda la lactancia.

  • Falsos positivos (destinar recursos extra a vacas que no rendirán) elevan los costos de alimentación y manejo sin retorno.

Nuestro Random Forest, al incorporar los dos primeros meses de litros, grasa, proteína y RCS, reduce estos errores donde importa de verdad: en este 25 % de alta variabilidad. De este modo:

  • Minimiza costos de oportunidad al no dejar escapar vacas valiosas.

  • Evita sobrecostos de manejo en animales de baja producción oculta.

  • Optimiza la rentabilidad del rodeo, concentrando recursos en los casos donde la simple regla de un mes no basta.

En definitiva, focalizar el modelo en este subgrupo de alta variabilidad se traduce en mejoras económicas reales y en una gestión más eficiente del tambo.

# Boxplot de la desviación estándar de producción mensual
ggplot(varianza, aes(y = sd_meses)) +
  geom_boxplot(fill = "#87CEFA", outlier.color = "red") +
  labs(
    title = "Boxplot de la desviación estándar de producción mensual",
    y     = "Desvío estándar (litros)",
    x     = ""
  ) +
  theme_minimal()

9. Reflexión final:

1. ¿Qué mido con la desviación estándar de producción mensual?

Al calcular el desvío estándar (sd_meses) de los primeros 12 meses de lactancia, obtuve:

  • Una mediana cercana a 8 L, lo que significa que la mitad de las vacas varía en torno a ±8 L mes a mes respecto a su propio promedio.
  • Un IQR de aproximadamente 5 – 10.6 L, es decir, el 50 % central de las curvas de producción permanece dentro de ese rango de dispersión.

Cuando observo vacas con un desvío muy bajo (por ejemplo, < 5 L), sé que su curva de producción es prácticamente plana. En esos casos, un único “mes piloto” basta para distinguir si serán de alta o baja producción.

2. ¿Qué me dice el umbral del 75 percentil?

Calculo el percentil 75 de sd_meses y obtengo ≈ 10.63 L. Eso implica:

  • El 25 % de mi población (unos 70 de 279 animales) presenta sd_meses > 10.63 L.
  • Estas vacas oscilan más de 10 L mes a mes, por lo que aplicar un único corte (“si en abril superó X L → alta”) resultaría muy impreciso: sus siguientes meses pueden subir o bajar demasiado.

3. ¿Por qué aporta valor el Random Forest?

  • Vacas de baja variabilidad (sd_meses ≤ 10.6 L; 75 % del total):

He comprobado que una regla univariable (por ejemplo, “litros abril > umbral”) predice casi con la misma exactitud que mi modelo complejo. En este subgrupo, el RF no mejora sustancialmente, pero tampoco perjudica.

  • Vacas de alta variabilidad (sd_meses > 10.6 L; 25 % restante):

No es sencillo predecir su producción usando un único mes de datos, pues su curva presenta fuertes irregularidades.. El Random Forest, en cambio, integra la información conjunta de los dos primeros meses (volumen, grasa, proteína y RCS) y logra reducir notablemente los falsos positivos y negativos comparado con cualquier criterio univariable.

4. Conclusión fundamentada

  • He comprobado que la variabilidad mes a mes marca el momento de usar un modelo multivariable.*
  • Si la curva es casi plana (desvío bajo), basta un mes piloto.
  • Si la curva presenta grandes oscilaciones (desvío alto), el Random Forest aprovecha múltiples señales tempranas para ofrecer una predicción mucho más robusta.”

De este modo, respaldo mi modelo con un criterio cuantitativo (la desviación estándar) y demuestro que su verdadero valor añadido se da justo en aquellas vacas cuya producción es suficientemente variable durante los primeros meses.

10.Conclusiones

En este trabajo diseñé y validé un modelo Random Forest capaz de anticipar, con datos de los dos primeros meses de lactancia y registros históricos, si una vaca alcanzará producción alta o baja a 305 días. El modelo supera ampliamente la línea base de predecir siempre la clase mayoritaria (Accuracy ≈ 89.2 % vs. 72.3 %, p-value ≈ 0.00017), con un Kappa de 0.73 que confirma un acuerdo sólido más allá del azar, alta sensibilidad para capturar vacas de alto rendimiento (93.3 %) y buena especificidad para identificar las de baja producción (78.3 %).

El análisis de importancia de variables revela que los volúmenes de leche de los dos primeros meses (lts_1, lts_2) y la producción histórica (leche_vital, leche_acum) son los predictores clave. Los parámetros de composición (grasa y proteína) aportan señal complementaria, mientras que el RCS y algunas variaciones mes a mes ejercen un rol menor. Esto sugiere enfocar la calidad de datos y el esfuerzo de recuento en los litros tempranos y, de manera secundaria, mantener aquellas variables de composición que realmente agreguen valor.

Para fundamentar por qué el modelo multivariable es especialmente útil, medimos la variabilidad mensual mediante la desviación estándar de producción (sd_meses). Solo un 25 % de las vacas presenta oscilaciones mayores a 10,6 L mes a mes; en ese subgrupo, una regla univariable basado en un solo mes sería imprecisa y generaría falsos positivos y negativos costosos. El Random Forest, al integrar múltiples mediciones tempranas, consigue corregir esos errores justo donde la curva de producción es más irregular, optimizando la asignación de recursos y minimizando pérdidas de litros o gastos innecesarios.

En síntesis, este modelo ofrece una herramienta práctica para tomar decisiones de manejo temprano:

  • Para la mayoría de vacas con baja variabilidad, una simple regla de umbral mensual ya predice con buena exactitud.

  • Para el 25 % con alta variabilidad, el Random Forest agrega verdadero valor al combinar lactancias tempranas, composición de leche y salud mamaria.

Recomiendo su aplicación como apoyo al tambero, complementado con monitoreos regulares y reentrenamientos periódicos, para maximizar la eficiencia y rentabilidad del rodeo.

11. Recursos adicionales

En el siguiente link se encuentra los dataset utilizados para el análisis, y el script para ejecutar en R.

Link: https://drive.google.com/drive/folders/1GJIltuAXF1TC_FRfRtbbvIZdN-wqMvkt?usp=sharing

.